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

1. Aug 16, 2015

### 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?

Code (C):
#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;
}

2. Aug 16, 2015

### Staff: Mentor

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= ...

Last edited: Aug 16, 2015
3. Aug 16, 2015

### Robin04

These are the equation I used.

#### Attached Files:

• ###### Spring_Pendulum_Sim.docx
File size:
47.9 KB
Views:
52
4. Aug 16, 2015

### FactChecker

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

5. Aug 16, 2015

### Robin04

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.

6. Aug 16, 2015

### 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?

7. Aug 16, 2015

### Robin04

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.

How do you mean damping?

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

8. Aug 16, 2015

### D H

Staff Emeritus
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.

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.

Last edited: Aug 16, 2015
9. Aug 16, 2015

### FactChecker

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.
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.)

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 .

10. Aug 16, 2015

### Robin04

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

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.

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

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?

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

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?

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

11. Aug 16, 2015

### Robin04

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

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

#### Attached Files:

• ###### Spring_Pendulum_Sim.docx
File size:
47.9 KB
Views:
42
12. Aug 16, 2015

### 256bits

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

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.

13. Aug 16, 2015

### 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 t2 + v0 t + y0.
At a time of t=0, that should be, with a = 9.81, a y-position of 4.91.

Does it?

14. Aug 16, 2015

### Robin04

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

15. Aug 16, 2015

### 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

16. Aug 16, 2015

### Robin04

I tried it, matched the data with my own calculations and they didn't agree. I realised 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

17. Aug 17, 2015

### D H

Staff Emeritus
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:
Code (C):
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
Code (C):
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
Code (C):

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:
Code (C):

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:
Code (C):
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:
Code (C):
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:
Code (C):

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

18. Aug 18, 2015

### 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.