Homework #8

Completing this assignment

For the each of the problems below, write a separate C++ program named in the problem problemN.cpp, where N is the problem number.

Important

Make sure your that g++ can compile your code. For some of the problems here, you will need to use -std=c++20.

Upload your C++ source files (not the executables produced by g++) to Brightspace.

Important

All work must be your own.

You may not use generative AI / large-language models for any part of this assignment.

1. File I/O : Update your projectile motion code from the previous homework (Homework #7) to write the output to a file. Have your main function do the integration for both \(C=0\) and \(C=0.3\), writing each integration to a separate file.

   solution

   Listing 163 projectile-rk-file.cpp

   #include <iostream>
   #include <vector>
   #include <cmath>
   #include <numbers>
   #include <format>
   #include <fstream>
   
   const double d = 0.074;
   const double A = std::numbers::pi * d * d;
   const double m = 0.145;
   const double rho_air = 1.2;
   const double g = -9.81;
   
   struct OrbitState {
       double t;
       double x;
       double y;
       double vx;
       double vy;
   };
   
   OrbitState rhs(const OrbitState& state, double C) {
   
       OrbitState dodt{};
   
       // dx/dt = vx; dy/dt = vy
   
       dodt.x = state.vx;
       dodt.y = state.vy;
   
       // d(vx)/dt = - 1/2 C rho_air A |v| vx / m
       // d(vy)/dt = - 1/2 C rho_air A |v| vy / m - |g|
   
       double v = std::sqrt(state.vx * state.vx + state.vy * state.vy);
   
       dodt.vx = -0.5 * C * rho_air * A * v * state.vx / m;
       dodt.vy = -0.5 * C * rho_air * A * v * state.vy / m - std::abs(g);
   
       return dodt;
   
   }
   
   OrbitState update_state(const OrbitState& state,
                           const double dt,
                           const OrbitState& state_derivs) {
   
       OrbitState state_new{};
   
       state_new.t = state.t + dt;
       state_new.x = state.x + dt * state_derivs.x;
       state_new.y = state.y + dt * state_derivs.y;
       state_new.vx = state.vx + dt * state_derivs.vx;
       state_new.vy = state.vy + dt * state_derivs.vy;
   
       return state_new;
   }
   
   void write_history(const std::vector<OrbitState>& history, double C) {
   
       std::string filename = std::format("projectile-C{:.2f}.out", C);
       std::ofstream ofile(filename);
       for (const auto& p : history) {
           ofile
               << std::format("{:11.8f} {:11.8f} {:11.8f} {:11.8f} {:11.8f}\n",
                              p.t, p.x, p.y, p.vx, p.vy);
       }
   }
   
   std::vector<OrbitState> integrate(const OrbitState& state0,
                                     const double C,
                                     const double dt_in) {
   
       // how the history of the orbit
   
       std::vector<OrbitState> orbit_history{};
       orbit_history.push_back(state0);
   
       double dt = dt_in;
   
       OrbitState state = state0;
   
       // integration loop
       while (state.y > 0.0) {
   
           // get the derivatives
           auto state_derivs = rhs(state, C);
   
           // get the derivatives at the midpoint in time
           auto state_star = update_state(state, 0.5 * dt, state_derivs);
           state_derivs = rhs(state_star, C);
   
           // update the state
           state = update_state(state, dt, state_derivs);
   
           orbit_history.push_back(state);
       }
   
       return orbit_history;
   
   }
   
   int main() {
   
       // set initial conditions
       OrbitState state0{};
   
       double C{};
       std::cout << "# Enter the drag coefficient C: ";
       std::cin >> C;
       std::cout << std::endl;
   
       // initial velocity and angle
       double v0 = 10.0;
       double theta = std::numbers::pi / 4.0;
   
       // initial height
       double y0 = 10.0;
   
       state0.t = 0.0;
       state0.x = 0.0;
       state0.y = y0;
       state0.vx = v0 * std::cos(theta);
       state0.vy = v0 * std::sin(theta);
   
       double dt = 0.05;
   
       auto orbit_history = integrate(state0, C, dt);
       write_history(orbit_history, C);
   
   }
   

   The only function that changed here is write_history, which now takes the drag coefficient, C, as an additional argument.

   In that function, the filename is constructed to include the drag coefficient, like:

   std::string filename = std::format("projectile-C{:.2f}.out"
   

   Note

   You could have also asked the user for the name of the filename or constructed it in some other fashion, as long as different values of C result in different files.

2. Lambda practice : Consider the following code:

   Listing 164 count_if.cpp

   #include <iostream>
   #include <algorithm>
   #include <vector>
   #include <cmath>
   
   bool is_perfect_square(int x) {
       int xroot = static_cast<int>(std::sqrt(x));
       return x == xroot * xroot;
   }
   
   int main() {
   
       std::vector<int> vec{1, 5, 9, 16, 21, 32, 36, 49, 60, 64};
   
       auto num = std::ranges::count_if(vec, is_perfect_square);
   
       std::cout << "there are " << num << " perfect squares in vec" << std::endl;
   
   }
   

   This counts how many elements in the vector vec are perfect squares, using the std::ranges::count_if algorithm.

   Rewrite this code to using a lambda function (review Lambda Functions) in the std::ranges::count_if call, in place of the is_perfect_square function above.

   solution

   Listing 165 count_if_lambda.cpp

   #include <iostream>
   #include <algorithm>
   #include <vector>
   #include <cmath>
   
   bool is_perfect_square(int x) {
       int xroot = static_cast<int>(std::sqrt(x));
       return x == xroot * xroot;
   }
   
   int main() {
   
       std::vector<int> vec{1, 5, 9, 16, 21, 32, 36, 49, 60, 64};
   
       auto num = std::ranges::count_if(vec, is_perfect_square);
   
       auto num2 = std::ranges::count_if(vec,
                                         [] (int x)
                                             {int xroot = static_cast<int>(std::sqrt(x));
                                              return x == xroot * xroot;});
   
       std::cout << "there are " << num << " perfect squares in vec" << std::endl;
       std::cout << "there are " << num2 << " perfect squares in vec" << std::endl;
   
   }
   

   In this case, our lambda function has two statements in the { }—that’s fine.

3. Transforming : std::ranges::transform takes a range (vector for us), an output iterator (where to start writing the result), and an operator (the function to apply to each element).

   For instance, given a vector v of double, we could do:

   std::ranges::transform(v, v.begin(), f)
   

   where f is a function of the form double f(double e), and the result would be to apply f(e) to each element, e, of v, updating our vector in-place.

   Let’s use this to convert a vector of indices into \(x\) values.

   Start with:

   std::vector<double> is{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
   

   and define \(x_\mathrm{min} = 1\), \(x_\mathrm{max} = 2\), and \(\Delta x = (x_\mathrm{max} - x_\mathrm{min}) / (N-1)\), where \(N\) is the number of elements in the vector.

   Now, apply the transformation: \(f(e) = x_\mathrm{min} + e \Delta x\) to the vector, using std::ranges::transform.

   Note

   You can use either a standard C++ function or a lambda-function for \(f(e)\).

   Finally, loop over the vector and output the updated elements.

   solution

   Listing 166 transform.cpp

   #include <algorithm>
   #include <iostream>
   #include <vector>
   
   int main() {
   
       std::vector<double> is{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
   
       double xmin{1.0};
       double xmax{2.0};
       double dx = (xmax - xmin) / (is.size() - 1);
   
       std::ranges::transform(is, is.begin(),
                              [=] (int i) {return i * dx + xmin;});
   
       for (auto e : is) {
           std::cout << e << std::endl;
       }
   
   }
   

4. Bounds : An operation we often want to do is search through a sorted list of numbers and find the interval that contains an desired value. This comes up in interpolation, for example.

   Consider the following vector:

   std::vector<double> temp_vec{0.1, 0.2, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 7.5, 10.0};
   

   For a value T0, we want to find the index into the vector, i, such that temp_vec[i-1] < T0 <= temp_vec[i]. We can get an iterator to temp_vec[i] using std::ranges::lower_bound.

   Note

   There is a similar function, std::ranges::upper_bound —they differ only in how they handle the case where we are searching for a value that is in our vector:

   • std::ranges::lower_bound returns an iterator to the first element \(\ge\) T0

   • std::ranges::upper_bound returns an iterator to the first element \(>\) T0

   Write a code that finds the index for:

   • T0 = 1.2

   • T0 = 3.0 —note: this is one of the data points in the vector.

   Tip

   You’ll want to use std::distance to convert the iterator into a distance from the beginning of the container, just like we did in class.

   Also check what happens in the case that our T0 is out of the limits of our temp_vec by finding the index for:

   • T0 = 0.05

   • T0 = 20

   solution

   Listing 167 binary_search.cpp

   #include <vector>
   #include <iostream>
   #include <algorithm>
   #include <format>
   
   int main() {
   
       std::vector<double> temp_vec{0.1, 0.2, 0.5, 1.0, 2.0,
                                    3.0, 4.0, 5.0, 7.5, 10.0};
   
       {
           // a point in-between elements in our vector
   
           double T0{1.2};
           auto it = std::ranges::lower_bound(temp_vec, T0);
           auto idx = std::distance(temp_vec.begin(), it);
   
           std::cout
               << std::format("{} is at index {}; in ({}, {}]\n",
                              T0, idx, temp_vec[idx-1], temp_vec[idx]);
       }
   
       {
           // a point matching an element in our vector
   
           double T0{3.0};
           auto it = std::ranges::lower_bound(temp_vec, T0);
           auto idx = std::distance(temp_vec.begin(), it);
   
           std::cout
               << std::format("{} is at index {}; in ({}, {}]\n",
                              T0, idx, temp_vec[idx-1], temp_vec[idx]);
       }
   
       {
           // a point less than any element in our vector
   
           double T0{0.05};
           auto it = std::ranges::lower_bound(temp_vec, T0);
           auto idx = std::distance(temp_vec.begin(), it);
   
           std::cout
               << std::format("{} is at index {}\n", T0, idx);
       }
   
       {
           // a point greater than any element in our vector
   
           double T0{20.0};
           auto it = std::ranges::lower_bound(temp_vec, T0);
           auto idx = std::distance(temp_vec.begin(), it);
   
           std::cout
               << std::format("{} is at index {}\n", T0, idx);
       }
   
   }
   

   Notice that for the case of 0.05, we get the index 0. So when we seek a value smaller than any element in our vector, we get the start of our vector.

   And for the case of 20.0, we get an index of 10 which is out of bounds for our vector. So when we seek a value larger than any element in our vector the iterator points to end().