Solving an ODE numerically with forward Euler method

The goal is to solve an ODE numerically with forward Euler method. The programs works well (numerical solution really near analytical one). The problem I see is that the Euler scheme don't jump to eyes, probably because of push_back() functions. This approach is the only one I found to let the time of simulation (number of steps) be variable (change only the t_max constant). Do you have any idea to improve the clarity of the program?

/*  R. M.
20.08.2012
Exercice 1.2 of Computational Physics, N. Giordano and H. Nakanishi
Euler method to solve: dv/dt = a - bv
*/

#include <iostream>
#include <fstream>
#include <string>
#include <vector>

void calculate(std::vector<double>& time, std::vector<double>& velocity, const double t_max, const double dt, const double a, const double b)
{
const int iterations(t_max / dt);

for(int i(1); i < iterations; i++)
{
time.push_back( time[i-1] + dt );
velocity.push_back( a * dt + (1 - b*dt) * velocity[i-1] );
}
}

void save(const std::vector<double>& time, const std::vector<double>& velocity, const std::string& filename)
{
std::ofstream file_out(filename);

for(int i(0); i < time.size(); i++)
{
file_out << time[i] << ' ' << velocity[i] << std::endl;
}

file_out.close();
}

int main()
{
// Constants
constexpr double t_max(8);          // Time to end simulation
constexpr double dt(0.01);          // Time step
constexpr double a(10.);            // Acceleration
constexpr double b(1.);             // Air friction

std::vector<double> time({0});      // Initial time (t = 0)
std::vector<double> velocity({0});  // Initial velocity (v = 0)

calculate(time, velocity, t_max, dt, a, b);
save(time, velocity, "veocity.dat");

return 0;
}


The smallest change that might make it more readable is to resize the vectors first, so you can index the elements rather than calling push_back.

It's generally a good idea to reserve vectors before a loop anyway if you know how big they're going to end up, and resize will handle that too.

void calculate(std::vector<double>& time, std::vector<double>& velocity,
const double t_max, const double dt, const double a, const double b)
{
const int iterations(t_max / dt);
time.resize(iterations, 0.);
velocity.resize(iterations, 0.);

for(int i(1); i < iterations; i++)
{
time[i] = time[i-1] + dt;
velocity[i] = a * dt + (1 - b*dt) * velocity[i-1];
}
}


if you're keeping two vectors synchronized like this, it's sometimes nicer to replace them with a single vector whose elements have two fields:

struct step {
double time;
double velocity;
};

void calculate(std::vector<step>& steps,
const double t_max, const double dt, const double a, const double b)
{
const int iterations(t_max / dt);
steps.resize(iterations);

for(int i(1); i < iterations; i++)
{
steps[i].time = steps[i-1].time + dt;
steps[i].velocity = a*dt + (1 - b*dt) * steps[i-1].velocity;
}
}


this is more of an aesthetic preference though, and it isn't clearly better in your case.

• I like very much your solution. Euler scheme is easier to read! Thank you! – user15816 Aug 21 '12 at 10:47