# Why Is My Runge Kutta Implementation Not Working for Car Simulation?

• platipo
In summary: C*sin(PI*OMEGA*t-t*τ), where τ is the time step. What you are doing is computing the acceleration of the lead car at every time step, which is u_n(t) = u_1(t)+C*sin(PI*OMEGA*t-t*τ). This is not a problem, as long as u_n is always computed, but it is not. Instead, u_n is computed only once, at the start of the simulation, and stored in global variable u.Issue #2You are using a Runge Kutta second order Method, but you want to make it a Runge Kutta
platipo
I've been browsing the forum, but I haven't found a way to solve this thing, which has been frustrating me for quite a while, I'm trying to develop a simulation of cue formation and propagation on cars on a single lane road, and of course, I have to compute the acceleration of vehicles obeying a certain car-following dynamic. I have a working Euler method, but I wanted to make it a Runge Kutta second order. And I can't... it's been two weeks now, I'm growing frustrated.
anyway, here's the "relevant" code, the function that I'm iterating over time.

Code:
void next(void){
int n;
double xw[CAR_NUMBER+1],uw[CAR_NUMBER+1];

double C=AMPLI/OMEGA;
u[CAR_NUMBER] = U0+AMPLI*cos(PI*OMEGA*t);
x[CAR_NUMBER] = U0*t+C*sin(PI*OMEGA*t);
// the first car in the cue is accelerating and braking periodically to perturbate the cue
Force(x,u,t+dt/2);
//the Force function returns a[n] which is the acceleration, I know it works, it gives no problem at all with the euler method
fprintf(fpx,"%10.3lf   ",a[CAR_NUMBER-1]);
for( n=0; n<CAR_NUMBER; n++ ){
xw[n]= x[n] +u[n]*dt/2;

uw[n]= u[n] +a[n]*dt/2;

}
//end of the trial step

Force(xw,uw,t+dt);
for( n=0; n<CAR_NUMBER; n++ ){
x[n]= x[n] +uw[n]*dt;
u[n]= u[n] + a[n]*dt;
}
//end of the second and final step
t=t+dt;
}

whenever I run this, the car after the first (which is of course just oscillating) brakes pretty hard, and quickly becomes stationary, with all the ones following reasonably doing the same. I Thank you in advance for any possible help, I found out about this forum just a couple of days ago and I find it really great. Sorry for my poor english, I'm doing my best...

You claim you know the function Force works. We don't. What is the time argument passed the Force, and why is this t+dt/2 for the intermediate step, t+dt for the final step?

sure; actually the force function now has no use for the time variable, but in a previous version the leader car motion was determined by the force function and that's why I'm passing time to that, it's actually a mistake, but I'm pretty sure it's not that relevant, in the program now, and anyway, as I had already said, I have a prefectly working euler method simulation using that force function. oh and btw (is btw tollerated as an abbreviation in this forum?) yes, t+dt is the final step and t+dt/2 is the half step.

Code:
void Force(double x_[], double u_[],double tau)
{

double d_r,d_s,uu,eps,K;
eps=0.3;
for( int n=0; n<CAR_NUMBER; n++ ) {
d_r= x_[n+1]-x_[n];
d_s= distanza(u_[n]);
if( d_r > d_s ) a[n]= -alfa*( u_[n]-(1+eps)*u_[n+1] );
else          a[n]= -alfa*( d_s- d_r );

}
}

then I guess you'll want to know what the "Distanza" function is doing, which is just retuning speed plus the minimum distance between cars, this was a first try, so I just made it linear.

thank you very much

The mid point method is essentially the same as dividing the time step by 2.

Since computation time isn't significant here, an implicit (predictor-corrector) iterative algorithm can be used. I explained a simple implicit algorithm here, which also includes a link to an Excel spreadsheet that does an 8 step iteration.

https://www.physicsforums.com/showpost.php?p=2285015&postcount=10

also here, showing the initial Euler step, but you can just start off with yg+0 = yn, and the first iterative step will be Euler.

http://en.wikipedia.org/wiki/Predictor-corrector_method#Euler_Trapezoidal_Example

Last edited:
platipo, your basic setup looks OK. I do have a comment about your code: You are doing a lot of communicating via global variables. That, in general, is not a particularly good idea.

Regarding your specific problem, here are two big issues:

Issue #1
platipo said:
Code:
     double C=AMPLI/OMEGA;
u[CAR_NUMBER] = U0+AMPLI*cos(PI*OMEGA*t);
x[CAR_NUMBER] = U0*t+C*sin(PI*OMEGA*t);
The lead car's position as a function of time is $x_n(t) = u_0 t + c\sin(\pi\omega t)$. Differentiating with respect to time, $u_n(t) = u_0 + c\pi\omega\cos(\pi\omega t)$. You have $u_n(t) = u_0 + c\omega\cos(\pi\omega t)$. Your lead car is violating the laws of physics. Moving to a smaller time step (the trial step in RK2) is going to make your simulation more sensitive to the error in your velocity calculation.

For that matter, why do you have this factor of PI?Issue #2
platipo said:
then I guess you'll want to know what the "Distanza" function is doing, which is just retuning speed plus the minimum distance between cars, this was a first try, so I just made it linear.
I'm going to exaggerate the problem here so you can see it. Rhetorical question: Does adding ten kilograms to fifteen meters make any sense? Does comparing ten kilograms to fifteen meters? Answer: Of course not. Length and mass are incommensurable quantities. Similarly, it doesn't make sense to add or compare a distance and a velocity. You should be adding and comparing commensurable quantities only. This is going to play havoc with your simulation when you reduce the time step (i.e., the first half step in RK2).
Now for a couple of questions.

Question 1: You mentioned your Euler propagation is working fine. How exactly did you compute that? In particular, which one of the following schemes did you use:

a. Update position, then velocity
Code:
     Force(xw,uw,t+dt);
for( n=0; n<CAR_NUMBER; n++ ){
x[n]= x[n] +u[n]*dt;
u[n]= u[n] + a[n]*dt;
}
b. Update velocity, then position
Code:
     Force(xw,uw,t+dt);
for( n=0; n<CAR_NUMBER; n++ ){
u[n]= u[n] + a[n]*dt;
x[n]= x[n] +u[n]*dt;
}

The first approach is *the* Euler method. The second is variously called semi-implicit Euler, Euler-Cromer, or symplectic Euler. That simple matter of switching which of position versus velocity is updated first makes an amazing huge difference in stability. The Euler method is not stable; symplectic Euler is.Question 2: What is your step size compared to the period 2/OMEGA? The period should be considerably larger (10 to 100 times larger) than the step size. In fact, you have what is called a stiff system. The step size should be tiny compared to the period.

D H said:
platipo, your basic setup looks OK. I do have a comment about your code: You are doing a lot of communicating via global variables. That, in general, is not a particularly good idea.
I know... I'm trying to improve my programming habits, bur sometimes I just get lazy.

D H said:
Issue #1
The lead car's position as a function of time is $x_n(t) = u_0 t + c\sin(\pi\omega t)$. Differentiating with respect to time, $u_n(t) = u_0 + c\pi\omega\cos(\pi\omega t)$. You have $u_n(t) = u_0 + c\omega\cos(\pi\omega t)$. Your lead car is violating the laws of physics. Moving to a smaller time step (the trial step in RK2) is going to make your simulation more sensitive to the error in your velocity calculation.

For that matter, why do you have this factor of PI?
no, C is ampli/omega, so the algebra was right, though as a matter of fact I have no use for that PI factor now, I'm removing it

D H said:
Issue #2
I'm going to exaggerate the problem here so you can see it. Rhetorical question: Does adding ten kilograms to fifteen meters make any sense? Does comparing ten kilograms to fifteen meters? Answer: Of course not. Length and mass are incommensurable quantities. Similarly, it doesn't make sense to add or compare a distance and a velocity. You should be adding and comparing commensurable quantities only. This is going to play havoc with your simulation when you reduce the time step (i.e., the first half step in RK2).
well, that doesn't make much sense, I'm not comparing quantities, I'm making the safe-distance a linear function of speed, and I'm adding to that a minimum distance

D H said:
Question 1: You mentioned your Euler propagation is working fine. How exactly did you compute that? In particular, which one of the following schemes did you use:

a. Update position, then velocity
Code:
     Force(xw,uw,t+dt);
for( n=0; n<CAR_NUMBER; n++ ){
x[n]= x[n] +u[n]*dt;
u[n]= u[n] + a[n]*dt;
}
b. Update velocity, then position
Code:
     Force(xw,uw,t+dt);
for( n=0; n<CAR_NUMBER; n++ ){
u[n]= u[n] + a[n]*dt;
x[n]= x[n] +u[n]*dt;
}

I used "the" Euler method, I'm trying the Symplectic now.

D H said:
Question 2: What is your step size compared to the period 2/OMEGA? The period should be considerably larger (10 to 100 times larger) than the step size. In fact, you have what is called a stiff system. The step size should be tiny compared to the period.

My step size is in the worst case scenario about 1/100 of the period, and I seriously hope that's tiny enough
now... why is my system still collapsing? any clue about what else could be going wrong?

Last edited:
platipo said:
D H said:
no, C is ampli/omega, so the algebra was right
No, the algebra is not right. I showed in my previous post that it isn't.

platipo said:
though as a matter of fact I have no use for that PI factor now, I'm removing it
Good. Does your system still have problems?

platipo said:
D H said:
Similarly, it doesn't make sense to add or compare a distance and a velocity. You should be adding and comparing commensurable quantities only.
well, that doesn't make much sense, I'm not comparing quantities, I'm making the safe-distance a linear function of speed, and I'm adding to that a minimum distance
You cannot add distance to speed. Doing so will create nonsense. Trying to compute the sum of 20 km/hr and 10 meters is just as meaningless as 10 meters + 15 kilograms. Make your diztanza function return a distance, not some meaningless nonsense. Until you fix this problem, all bets are off regarding the performance of your program.

Once you fix this basic problem (and you really need to fix this problem), you should ask yourself whether this algorithm reflects reality. Think of it this way. Suppose you are driving just behind the lead car and that a long, long line of cars is behind you. Suppose the 10000th and 10001st cars in the just got a little bit too close together. Realistically, will that make you slow down? You are probably much better off to consider just the distances between your car and the car immediately in front of you and between your car and the car immediately behind you.

platipo said:
now... why is my system still collapsing? any clue about what else could be going wrong?
I already told you but you refused to listen. Fix the velocity, and fix your distanza function. Make it return a distance, not a jumble of distance and speed.

D H said:
Once you fix this basic problem (and you really need to fix this problem), you should ask yourself whether this algorithm reflects reality. Think of it this way. Suppose you are driving just behind the lead car and that a long, long line of cars is behind you. Suppose the 10000th and 10001st cars in the just got a little bit too close together. Realistically, will that make you slow down? You are probably much better off to consider just the distances between your car and the car immediately in front of you and between your car and the car immediately behind you.

numbers in a calculator are just numbers; I want the safe distance to be a linear function of speed, and that's what I'm getting. How would you do that?
The algorithm DOES reflect reality (somehow) it's called a "car following" model, and there's plenty of literature about that, and you're right, in reality you would bother only about the car preceding you and the one following you, in this simple model you are bothering only about the car preceding you, if I'm just behind the leader in this program you care only about the changes in distance and relative velocity with that car, that's what the program is doing.
my program is working way better now anyway, I've changed a few things, I have of course some problems with the stiffness, but I'm trying to work my way trough that... I'll keep this thread updated, thank you very much for your help.

platipo said:
numbers in a calculator are just numbers; I want the safe distance to be a linear function of speed, and that's what I'm getting.
That is not what you are getting. You are getting nonsense.

How would you do that?
I might use something like the rule taught in drivers education. What is the safe following distance rule?

D H said:
That is not what you are getting. You are getting nonsense.

I might use something like the rule taught in drivers education. What is the safe following distance rule?

they say the proper formula is this
$$s = v^2 / (2 · g · k)$$
BUT I didn't want to be a square function of speed but a linear one, and putting all constants equal to one... this is not nonsense of course, everything might make sense dimensionally provided you multiply it by one the right number of times right?

Try using the three second rule.

D H said:
Try using the three second rule.
well, actually we have many different safe distance models we are working on, and the dynamics itself is all but done, there are many different models that are worth a try, the one I'm using is one, it doesn't take into account reaction times and it's linear, but before switching to a more complex model I wanted the model itself to react deterministically using a reasonable ODE solving method, and the euler method isn't really an option, now the Runge Kutta is working, so I'll see... thank you again.

I have it! It finally works! Thank you again, I'll soon try out some different dynamics and see what I can get out of it.

## 1. What is a Runge-Kutta method?

A Runge-Kutta method is a numerical method used to solve ordinary differential equations (ODEs). It is an iterative algorithm that approximates the solution of an ODE by breaking it down into smaller steps.

## 2. Why is this thread called "Another dumb Runge Kutta thread"?

The title of this thread is meant to be humorous and playful. It is not meant to disparage the Runge-Kutta method, but rather to poke fun at the fact that there have been many discussions and debates about this topic.

## 3. Is the Runge-Kutta method the most accurate method for solving ODEs?

No, there are other numerical methods such as the Euler method and the Adams-Bashforth method that may be more accurate in certain situations. The choice of method depends on the specific ODE and the desired level of accuracy.

## 4. What are the limitations of the Runge-Kutta method?

One limitation of the Runge-Kutta method is that it can be computationally expensive for complex ODEs. It also may not be suitable for stiff equations, which are ODEs that involve rapidly changing values.

## 5. Can the Runge-Kutta method be applied to systems of differential equations?

Yes, the Runge-Kutta method can be extended to solve systems of differential equations. It involves using the same principles but with multiple equations and variables.

• Programming and Computer Science
Replies
15
Views
2K
• Differential Equations
Replies
65
Views
3K
• Programming and Computer Science
Replies
1
Views
2K
• Programming and Computer Science
Replies
5
Views
2K
• Programming and Computer Science
Replies
14
Views
5K
• Programming and Computer Science
Replies
2
Views
3K
• Programming and Computer Science
Replies
4
Views
7K
• Mechanics
Replies
9
Views
2K
• Differential Equations
Replies
6
Views
3K
• Programming and Computer Science
Replies
8
Views
13K