3
\$\begingroup\$

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;
}
\$\endgroup\$
2
\$\begingroup\$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy