Using an external array in ode45

I am trying to use an external array for use in ode45. The code should be like for the first iteration of the ode45, it picks up the first element and for the second iteration it picks up the second one. I have limited the size of ode45 using a predetermined time array (right now upto 200 iterations for 10 seconds). Can someone guide me how to do this?

I tried using an external variable "n" which was set to value 1. It would get incremented upon the finishing of the iteration but on the next iteration, it would be set to 1 again. I even used a for loop and got the same results.

[Jan]

The question is not clear.

"I have limited the size of ode45 using a predetermined time array" - do you mean the input tspan? This determines the number of output points, not the number of steps or of calls of the function to be integrated.

What do you call "first iteration"? What is iterated? What do you call "an external variable 'n'"?

Please post your code.

[Kashish Pilyal]

Yes, the input tspan. By "first iteration", I mean the first calculation the ode45 solver does. It calculates the values in the function and integrates them overtime. Then it uses those values as input in the next calculation. I hope I am correct about how ode45 works.

This is the main code

%% Leader vehicle
q=0; %position
v=-2; %velocity
a=0; %acceleration
x_0=[q;
   v;
   a]; %initial conditions
t=linspace(0,30,20);    
L1=4; %lenght of car(m)
[t,y]=ode45(@leader,t,x_0);   

%% for the first vehicle (leader)
figure(1)
plot(t,y(:,1));
hold on
plot(t,y(:,2));
hold on
plot(t,y(:,3));

%% for the followers

dr=5; %gap b/w cars(m)

%first follower, second vehicle
q1=L1+dr;
v1=0;
a1=0;


n=1;
x1=[q1;
    v1;
    a1];
[t,y1]=ode45(@(t,x)follower(t,x,y,n),t,x1);

The leader code

function xdot= leader(t,x)
%% Parameters

tau_i=0.1;
h_i=0.5;


k_d=0.7;
k_p=0.2;
T=tau_i/h_i;



%% u_i and eta and errors

%for leader vehicle
e_i1=0;
e_i2=0;
   

u_i=T*[k_p k_d]*[e_i1;e_i2]+ ( (1-T)*x(3) );

%% CACC equations
x_1dot=x(2); %q(position)dot
x_2dot=x(3); %v(velocity)dot
x_3dot=(-(1/tau_i)*x(3))+((1/tau_i)*u_i); %a(accelration)dot

xdot=[x_1dot;
      x_2dot;
      x_3dot];

and the follower code

function xdot=follower(t,x,x_p,n)

tau_i=0.1;
h_i=0.5;


k_d=0.7;
k_p=0.2;
T=tau_i/h_i;


q_p=x_p(n,1);
v_p=x_p(n,2);
a_p=x_p(n,3);

%%
e_i1=q_p-(  x(1)+(h_i*x(2))  );
e_i2=v_p-(  x(2)+(h_i*x(3))  );
    
u_i=T*[k_p k_d]*[e_i1;e_i2]+ ( (1-T)*x(3) )+ ( T*a_p );
    
    
%% CACC equations
x_1dot=x(2); %q(position)dot
x_2dot=x(3); %v(velocity)dot
x_3dot=(-(1/tau_i)*x(3))+((1/tau_i)*u_i); %a(accelration)dot

xdot=[x_1dot;
      x_2dot;
      x_3dot];
 
 n=n+1;

As you can see that I have used a variable "n" to externally give a value and increase its value according to the iteration of the calculation it does.

Answer 1

[Jan]

This cannot work for several reasons:

1. ODE45 is designe to integrate smooth functions with a certain accuracy. Using non-smooth parameters let the function follower() be non-differentiable. Therefore the stepsize controller of ODE45 can drive mad. You might get a final value, but there is a chance, that it is dominated by rounding errors.
2. ODE45 rejects steps, if the error tolerance is not met. Then a step is repeated with a smaller step size. If this happens, your n is still increased.
3. The number of elements in tspan does not mean, that the function to be integrated is called the same number of times. For one step of the integration, the function to be integrated is called multiple times. The step size is determined automatically controlled by the given tolerances. So the output it the steps of tspan should be created by an interpolation of the trajectory.
4. n is provided as input argument. Then increaing n inside the function does not change the value outside the function magically. If this would be the case, programming would force the user to know all internally used variables of a function to avoid collisions.

You need a completely different design.

Answer 2

[Walter Roberson]

The mathematics of the Runge-Kutta methods require that the function and two derivatives of the function are continuous. When you modify the parameters of the function on the fly, you make the function discontinuous. If you are lucky ode45 will notice and quit the function. If you are not lucky, ode45 will just return incorrect results.

Also, as Jan indicated, ode45() calls the function multiple times at each location. Typically it calls 6 times per at each location and proposed step size. If it decides that the change in function values is too large, it will reject the proposal and try again with a smaller step. No information is passed to the ode function about which sub-evaluation is being done, or about whether a step was rejected or not. ode45 and some of the other RK functions sometimes go backwards in time even without having rejected anything. (Also, the points that are actually output are generally interpolated: when it notices that it has successfully passed one of the times output is requested for, then it will interpolate the results for the requested time.)

So... Your approach is all wrong.

In your situation what might possibly work is to use interp1 to interpolate the value of the variable over time. However you need an interpolation of at least degree 3, such as spline: the default linear interpolation will not satisfy the continuity requirements.

[Walter Roberson]

Also note that if your variable represents discrete events such as impulses or injection of chemicals, then you cannot use interpolation.

[Kashish Pilyal]

Thank you for the help. I am going to try another approach in which I will make the leader and follower code into one function so that it is continous and not changed on the fly. I will run that function as ode45. Will that work in my situation?