[C++] What's wrong with my program?

[Robin04]

Hi,

I made a small program that simulates the motion of a body attached to an ideal spring and released from a horizontal position. Some more info here (especially #1 and #3): https://www.physicsforums.com/threads/equations-of-motion-unsolvable-with-elementary-method.827826/

I've been thinking a lot, I've already checked several times my equations and my code but I don't know what's wrong. The body just flies away to the other side of the Sun. :D

Do you have any idea what's the problem with my code?

#include <iostream>
#include <math.h>
#include <fstream>
using namespace std;

int main()
{
    float    x = 0.0f, y = 0.0f,
            v_x = 0.0f, v_y = 0.0f,
            a_x = 0.0f, a_y = 0.0f,
            D = 0.0f, L = 0.0f, m = 0.0f,
            dt = 0.0f, T = 0.0f, t = 0.0f, alpha = 0.0f, F = 0.0f;

    ofstream File;

    cout << "D = ";
    cin >> D;
    cout << "L = ";
    cin >> L;
    cout << "m = ";
    cin >> m;
    cout << "dt = ";
    cin >> dt;
    cout << "T = ";
    cin >> T;

    File.open("Coordinates.txt", ios::out);
    File << "D = " << D << "\tL = " << L << "\tm = " << m << "\tdt = " << dt << "\tT = " << T << endl;
    x = L;  // The pendulum starts from a horizontal position; y = 0

    while (t <= T)
    {

        if (y != 0)
            alpha = atan2(x, y); //avoiding dividing by 0
        else
            alpha = 3.1415f / 2.0f; //if y = 0, than the angle is 90 degrees

        F = D*(hypotf(x, y) - L); //calculating the spring force: D * dl = D * (sqrt(x^2 + y^2) - L)

        a_x = F / m * sin(alpha); //calculating the components of the acceleration
        a_y = 9.81f - F / m * cos(alpha);

        cout << "Alpha = " << alpha*180/3.1415f << endl; //These are for trouble shooting to see what values go wrong, and as everything depends on everything this brought no help
        cout << "dl = " << hypot(x, y) << endl;
        cout << "F = " << F << endl;
        cout << "a_x = " << a_x << endl;
        cout << "a_y = " << a_y << endl;

        v_x += (a_x * dt); //velocities
        v_y += (a_y * dt);

        cout << "v_y = " << v_x << endl;

        x += (0.5f * a_x * pow(dt, 2) + v_x * dt); //updating the positions
        y += (0.5f * a_y * pow(dt, 2) + v_y * dt);

        cout << "t = " << t << "\tx = " << x << "\ty = " << y << endl;
        File << t << "\t" << x << "\t" << y << endl;

        t += dt;
    }

    File.close();

    system("pause");
    return 0;
}

 

[jedishrfu]

Have you worked out the problem on a sheet of paper with columns for x, y, vx, vy, ax, ay ...?

From there you can see where things go astray.

Perhaps the order of your equations matters like try:

new x= old ax and vx used
new y= old ay and vy used
new vx= old ax used
new vy= old ay used
new ax= ...
new ay= ...

 

[Robin04]

These are the equation I used.

 

[FactChecker]

dt = 0? Try dt = 0.01; // time step

 

[Robin04]

dt = 0? Try dt = 0.01; // time step

That's just an initial value. It's my habit to set every variable to 0 in order to avoid using them without value (in case I forgot to set them). The actual dt the program uses comes from an input at the beginning of the program.

 

[FactChecker]

Sorry, I didn't notice the initialization of dt.
1) It sounds like the sign of force is not correct if it flies off to infinity.
2) You might want to put some damping in your equations.
3) I see that you omitted gravity from your forces. Was that intentional?

 

[Robin04]

1) It sounds like the sign of force is not correct if it flies off to infinity.

Yes, I already thought about that and tried it with the opposite direction but it didn't work. By the way according to the equations in the doc, the signs are correct.

2) You might want to put some damping in your equations.

How do you mean damping?

3) I see that you omitted gravity from your forces. Was that intentional?

It's there in the a_y. It start's by 9.81f -...

 

[D H]

Aside: Yes, I am being mean and nasty. Please do not take it personally. You asked what's wrong with your code. I'm treating your code as if you had asked me to peer-review it. I'm mean and nasty in peer reviews. I don't pick at the person, but I definitely do pick at the code they wrote. I want exactly the same (actually, more so) in return. I want other people to be extremely mean and extremely nasty when it comes to reviewing the code that I myself have written.

        F = D*(hypotf(x, y) - L); //calculating the spring force: D * dl = D * (sqrt(x^2 + y^2) - L)

        a_x = F / m * sin(alpha); //calculating the components of the acceleration
        a_y = 9.81f - F / m * cos(alpha);

The above lines of code are the immediate source of your problem. You have gravity plus a spring in the y-component (which is fine), but you have an anti-spring in the x-component (which is anything but fine). You have a sign error in the calculation of a_x. The immediate fix is to make that x-coordinate anti-spring force into a spring force by negating the sign: a_x = -F / m * sin(alpha) rather than a_x = F / m * sin(alpha).

However, (or perhaps However ), That is not the end of your problems. Why do you have that unneeded protection against y == 0 in your calculation of alpha =atan2(x, y)? For that matter, why are you using this construct at all? For that matter, why are you using alpha, period? (It's not needed.) Why are you using float rather than double? And (this is a bit more advanced), why are you using symplectic Euler integration?Your basic problems will disappear if you fix that one line (change a_x = F / m * sin(alpha) to a_x = - F / m * sin(alpha).

With regard to your deeper problems: I'll get to that in a couple of hours or so. My dogs (see my avatar) haven't had their long, long walk for six days. The weather here has cooled from nasty, nasty hot&humid to just mildly hot & mildly humid. My puppies need their long, long walk, and so do I.

 

[FactChecker]

Yes, I already thought about that and tried it with the opposite direction but it didn't work. By the way according to the equations in the doc, the signs are correct.

Ok. Check the signs of force, signs of acceleration, change of velocity, rate of change of position, in that order. If it is not returning back, the first one that is wrong will tell you where the problem is.

How do you mean damping?

In motion through air or a fluid, the velocity creates an opposing force that is proportional to velocity squared. Anything that is not in outer space sees some resistance to motion (friction, air, fluid, etc.)

It's there in the a_y. It start's by 9.81f -...

Sorry. I missed that.

PS. I can't argue with D.H.s comments, but all-in-all your code is a lot better than a lot of code I see. The variable names make sense. There are comments (could be more obvious to divide sections. Like a top level outline.). It is logically organized .

 

[Robin04]

Aside: Yes, I am being mean and nasty. Please do not take it personally. You asked what's wrong with your code. I'm treating your code as if you had asked me to peer-review it. I'm mean and nasty in peer reviews. I don't pick at the person, but I definitely do pick at the code they wrote. I want exactly the same (actually, more so) in return. I want other people to be extremely mean and extremely nasty when it comes to reviewing the code that I myself have written.

That's not a problem at all. In fact, I completely agree with you. I can learn the most from straightforward opinions.

The above lines of code are the immediate source of your problem. You have gravity plus a spring in the y-component (which is fine), but you have an anti-spring in the x-component (which is anything but fine). You have a sign error in the calculation of a_x. The immediate fix is to make that x-coordinate anti-spring force into a spring force by negating the sign: a_x = -F / m * sin(alpha) rather than a_x = F / m * sin(alpha).

I tried it but it didn't work. The body still flies to space. Actually I don't really understand why I have a sign error. I checked my equations several times and they seem fine to me, but there has to be something that I miss. I uploaded the equations in #3.

Why do you have that unneeded protection against y == 0 in your calculation of alpha =atan2(x, y)

I accidentally left it there. In an earlier version I had problems with the y coordinates, the loop started from y = 0.

For that matter, why are you using alpha, period? (It's not needed.)

My reason was that it was simpler to just copy the equations to the code and thus it's easier to check for any errors. By the way how can I determine the components of the force without knowing the angle between the vector and an axis?

Why are you using float rather than double?

I haven't programmed for a long time and I forgot some things. I corrected them. :)

why are you using symplectic Euler integration?

As I'm a high school student I often don't know what I'm exactly doing, I just do what makes the most sense to me, or what I've been taught. But this fact doesn't mean that I don't even want to know what I"m doing. :D What method do you think would be the most accurate?

With regard to your deeper problems: I'll get to that in a couple of hours or so. My dogs (see my avatar) haven't had their long, long walk for six days. The weather here has cooled from nasty, nasty hot&humid to just mildly hot & mildly humid. My puppies need their long, long walk, and so do I.

Have fun for me too. Here we still have 38 degrees. :P :D

 

[Robin04]

These are the equation I used.

I accidentally left out a cos alpha in the y direction forces, but the accelerations were correct. I uploaded the new doc.

In motion through air or a fluid, the velocity creates an opposing force that is proportional to velocity squared. Anything that is not in outer space sees some resistance to motion (friction, air, fluid, etc.)

Yes, I thought about putting air resistance into the calculations but first I would like to solve this problem. :)

 

[256bits]

I haven't programmed for a long time and I forgot some things.

You can format your fields in the file to look nicer with setw() and setprecision(), and to make it easier to add or delete fields separate the output variables.

File << setw(width) << t;
File << setw(width+5) << setprecision(3) << alpha; // here is an added output variable
File << setw(width+5) << setprecision(3) << x; // one way to change the width by adding a constant value
File << setw(width+5) << setprecision(3) << y;
File << endl;

Add the line
const int width = 6; // field width - change the value as desired.
somewhere top of the code. Below the main(), or the initializtion would be OK.

 

[256bits]

Also, try setting D to 0 and see what your results are.
With no spring constant acting the mass should fall in only the y-direction.

Note that the values of computed y-position should come out following the equation y = 0.5 a t² + v₀ t + y₀.
At a time of t=0, that should be, with a = 9.81, a y-position of 4.91.

Does it?

 

[Robin04]

You can format your fields in the file to look nicer with setw() and setprecision(), and to make it easier to add or delete fields separate the output variables.

File << setw(width) << t;
File << setw(width+5) << setprecision(3) << alpha; // here is an added output variable
File << setw(width+5) << setprecision(3) << x; // one way to change the width by adding a constant value
File << setw(width+5) << setprecision(3) << y;
File << endl;

Add the line
const int width = 6; // field width - change the value as desired.
somewhere top of the code. Below the main(), or the initializtion would be OK.

Thank you very much. That's a very useful tool. :)

 

[256bits]

In addition, you can try the following:

File << setiosflags(ios::fixed);

File << setw(width) << t;
File << setw(width+2) << setprecision(3) << alpha;
File << setw(width+5) << setprecision(3) << x;
File << setw(width+5) << setprecision(3) << y;
File << endl;

Include file is,
#include <iomanip> // std::setiosflags

 

[Robin04]

Also, try setting D to 0 and see what your results are.
With no spring constant acting the mass should fall in only the y-direction.

Note that the values of computed y-position should come out following the equation y = 0.5 a t² + v₀ t + y₀.
At a time of t=0, that should be, with a = 9.81, a y-position of 4.91.

Does it?

I tried it, matched the data with my own calculations and they didn't agree. I realized that in this order in the first loop the program already calculates with initial velocities, which is not correct. So I rearranged the order of equations starting by calculating the positions, then the velocities and finally the accelerations and it worked, but only with D = 0. If I use this with D > 0 then the first loop will also be a free fall which is not the case (of course it's really close because at t=0 the spring is unloaded). But the flying to the Sun still happens. :D

 

[D H]

I tried it but it didn't work. The body still flies to space. Actually I don't really understand why I have a sign error. I checked my equations several times and they seem fine to me, but there has to be something that I miss. I uploaded the equations in #3.

Regarding your equations, I, for one, am not going to open a .docx file posted by an unknown person. Doing so is a recipe for a massive viral infection. Nothing against you personally, but I don't know you, so I can't trust you. This site uses MathJax to implement LaTeX. For example, $F = -D(r - L)$. Using LaTeX is the recommended practice at this site (and many others) for representing equations.

To see that you indeed do have a sign error, suppose the initial state is x=2L, y=0. Your code has the spring force away from rather than toward the center.

You don't need to find $\alpha$ to compute the acceleration:

r = hypot(x, y);
a = D/m*r*(1.0-L/r);
a_x = -x/r * a;
a_y = -y/r * a + 9.81;

Note that the multiplication and division by r can be factored out, resulting in the simpler form

r = hypot(x, y);
a_factor = D/m*(1.0-L/r);
a_x = -x * a_factor;
a_y = -y * a_factor + 9.81;

Another problem with your code is the manner in which you update the state (position and velocity). You are using

v_x += a_x * dt;
v_y += a_y * dt;
x += 0.5f * a_x * pow(dt, 2) + v_x * dt;
y += 0.5f * a_y * pow(dt, 2) + v_y * dt;

This is incorrect. I'll give a few simple examples of how to perform this numerical integration.

1. Basic Euler method.
A very simple way to integrate a first order ordinary differential equation, $\frac{dx}{dt} = f(x,t)$ is via Euler's method:

x += f(x,t)*dt;

Yours is a second order differential equation. Any nth order ordinary differential equation can be re-expressed as a first order ordinary differential equation by adding extra states. In this case, the extra states are the x and y components of velocity:
$$\frac {d}{dt} \left(\begin{matrix} x\\y\\v_x\\v_y\end{matrix}\right) =
\left(\begin{matrix} v_x\\v_y\\-D/m (1-\frac L r) x\\-D/m (1-\frac L r) y + 9.81\end{matrix}\right)$$

This leads to the following code:

while (not done) {
    // Calculate a_x and a_y as above
    x += v_x * dt;
    y += v_y * dt;
    v_x += a_x * dt;
    v_y += a_y * dt;
}

While you should never use the basic Euler method, you should understand how it works.1. Symplectic Euler method.
One problem with using the basic Euler method on a second order ODE is that this technique does not conserve energy. Simply switching the order in which position and velocity are updated comes close to conserving energy:

while (not done) {
    // Calculate a_x and a_y as above
    v_x += a_x * dt;
    v_y += a_y * dt;
    x += v_x * dt;
    y += v_y * dt;
}

I'm not going to go into the mathematics that makes this simple switch a vast improvement over the basic Euler method.

2. Leapfrog.
One problem with *any* Euler method is that they are rather inaccurate. The error is proportional to $\Delta t$. You need to decrease the step size to improve accuracy, but decreasing the step size runs into numerical representation issues. (For example, try adding 1.0 + 1e-16 on your computer. The result will be exactly 1.0). A slight change yields a technique that is accurate to second order (error is proportional to $\Delta t^2$), the leapfrog technique:

// Calculate accelerations based on initial position, e.g., a_x = 0.0, a_y = 9.81.
dv_x = 0.5*a_x * dt;
dv_y = 0.5*a_y * dt;
while (not done) {
    v_x += dv_x;
    v_y += dv_y;
    x += (v_x + dv_x * dt) * dt;
    y += (v_x + dv_y * dt) * dt;
    // Calculate acceleration as above, using the end-of-interval x and y
    dv_x = 0.5*a_x * dt;
    dv_y = 0.5*a_y * dt;
    v_x += dv_x;
    v_y += dv_y;
}

 

[Robin04]

Thank you very much for your answer. We have a power cut due to yesterday's storm. When the electricity comes back I'll try your suggested modifications and write back how they worked or if everything is clear.