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

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

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

 

[D H]

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?

Show your function Force, please.

 

[platipo]

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.

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

 

[rcgldr]

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 y_g+0 = y_n, and the first iterative step will be Euler.

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

 

[D H]

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

     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

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

     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

     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.

 

[platipo]

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.

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

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

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

     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

     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.

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
thank you very much for your time and your help.
now... why is my system still collapsing? any clue about what else could be going wrong?

 

[D H]

Your lead car is violating the laws of physics.

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.

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?

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 10000^th and 10001^st 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.

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.

 

[platipo]

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 10000^th and 10001^st 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.

 

[D H]

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?

 

[platipo]

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?

 

[D H]

Try using the three second rule.

 

[platipo]

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.

 

[platipo]

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.