Solving an ODE numerically with forward Euler method

Question

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;
}

Answer 1

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.