ode45 and rungekutta yield different result

I am using ode45 and rungekutta script to solve some piecewise differential equations. But the result is a bit different.

Can anyone help me to solve the problem? I am hoping to get a similar result as ode45.

Please refer to Torsten's reply for the answer.

Thanks,

[Sam Chak]

@haohaoxuexi1

Didn't check your dynamical equations due to code clutter issue.

Is it a piecewise smooth dynamical system, or a piecewise discontinuous system?

If your rungekutta code is 100% correct and is expected to return output similar to the ode45 solver, then I guess something is missing in the dynamical system.

[haohaoxuexi1]

it is a piecewise smooth system, the problem is the first peak looks correct, after the first peak, the value of the rungekutta one starts to go lower.

[haohaoxuexi1]

@Sam Chak would you mind to check it for me? I didn't see anything wrong in the code.

[Torsten]

@haohaoxuexi1

First change your Runge-Kutta code such that both ODE45 and your code can use the same function "SDOFFriction" to compute the derivatives.

Then someone might invest the effort to compare the codes.

[haohaoxuexi1]

@Torsten

I didn't get what you mean, can you please show me a simple example how to do this?

I am new to this.

Answer 1

[Torsten]

Both results look the same for me.

C_p=3.026700000000000e-10;
kp=392400000;
theta_p = 0.231516000000000;
R_s=20*1.8e6; % ohm
ks=1.6e3;
k=(kp*ks)/(kp+ks); 
m=0.1;
xi1=0.000001;
c=2*xi1*sqrt(k*m); 
v_band=2;
uD=0.2;
uS=0.5; 
C=20000; 
Fn=100;
% initial condition
x0=[0; v_band; 0];
dt=1e-4; % time step
tend=3; % time final
zeit=0:dt:tend; % 
Y = runge_kutta_RK4(@(t,y)SDOFFriction(t, y, m, c, k, v_band, uD, uS, C, Fn, theta_p, ks, kp, R_s, C_p),zeit,x0);
figure(1);
plot(zeit, Y(3,:),'r'); %voltage
box on;
grid on;
xlabel('Time [s]', 'fontsize', 12, 'FontWeight', 'bold', 'fontname', 'Times New Roman');
ylabel('Voltage [V]', 'fontsize', 12, 'FontWeight', 'bold', 'fontname', 'Times New Roman');

%%
options = odeset('RelTol',1e-9);
[T,Y] = ode45(@(t,y)SDOFFriction(t, y, m, c, k, v_band, uD, uS, C, Fn, theta_p, ks, kp, R_s, C_p),zeit,x0,options);
figure(2)
plot(T,Y(:,3))

function dx = SDOFFriction(t, x, m, c, k, v_band, uD, uS, C, Fn, theta_p, ks, kp, R_s, C_p) %Reibung

%t
vrel=v_band-x(2); % rel velocity
if abs(vrel)/v_band<0.001  %          
    fstatic = uS*Fn*sign(vrel); faply = k*x(1)+c*x(2)+theta_p*(ks/(kp+ks))*x(3);
    if abs(faply)<uS*Fn
        frictionf=k*x(1)+c*x(2)+theta_p*(ks/(kp+ks))*x(3);            
    elseif abs(faply)>=uS*Fn
        frictionf=fstatic;
    else                
    end
    dx(1)=x(2);
    dx(2)=(-c.*x(2)-k.*x(1)+frictionf-theta_p*(ks/(kp+ks))*x(3))/m;
    dx(3)=(theta_p*(ks/(kp+ks))*x(2)*R_s-x(3))/R_s/C_p;
else %
    mu=uD+(uS-uD)*exp(-(C*abs(vrel))); % 
    frictionf=mu*Fn*sign(vrel);
    dx(1)=x(2);
    dx(2)=(-c.*x(2)-k.*x(1)+frictionf-theta_p*(ks/(kp+ks))*x(3))/m;
    dx(3)=(theta_p*(ks/(kp+ks))*x(2)*R_s-x(3))/R_s/C_p;
end     

dx = dx'; % output result

end
%%

function Y = runge_kutta_RK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(n,N);
Y(:,1) = Y0;
for i = 2:N
    t = T(i-1);
    y = Y(:,i-1);
    h  = T(i) - T(i-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(:,i) = y + h/6*(k0+2*k1+2*k2+k3);
end
end 

[Torsten]

I gave you the correct and efficient form of a Runge-Kutta-code of order 4, and the results for your problem were the same as with ODE45.

If you want to use your code and find the error in it, you can do so. But don't you see that it's programmed completely inefficient ? Imaginge you had to solve 50 ODEs. Would you write a function handle for each of them ? I think you must admit that is nearly impossible to debug this ansatz already for your 4 ODEs.

One error in your code is the line

fyy1 = @(t,y1,yy1,Vr) (-c*yy1-k*y1+frictionf-(ks/(kp+ks))*theta_p*Vr)/m;

"frictionf" depends on all three unknown functions, and you had to include this expression here instead of the constant value "frictionf".

Further it must read

frictionf = frictionf_ary(i) 

instead of

frictionf = frictionf_ary(i+1)

[haohaoxuexi1]

@Torsten I am very thankful for you. And for sure, you code is simple and effecient as you said. It matched with ode45 because it used the same function file to do the calculation.

But I am asking if you can quickly go through my file for a comparison (the logic), I have spent quite a while on my file to find the error. What I can see if the first peak match, but after the first peak, it starts to misalign, and I want to know the reason.

[Torsten]

Look at my comment above and the errors I found.

It will be quite hard to repair the first one within your Runge-Kutta code form.

Such an error will always happen if you calculate expressions in your function that depend on the solution variables. All such expressions had to be included in your function handles in their written-out form.

Using one function file for all integrators you want to test in which - given t and the vector of solution variables y - you calculate and return the vector of derivatives dy/dt is the best way to circumvent any such kind of errors.

[Torsten]

I thought rungekutta is doing iteration calculation as well?

That's exactly the point. During one time step, Runge-Kutta evaluates your function for different values of y1, yy1 and Vr. But your value for "frictionf" (which also depends on y1, yy1 and Vr and should therefore change) is not changed accordingly because you include it as a constant in the function handle fyy1 with value frictionf_ary(i) at time t(i).

By the way: why did you delete your Runge-Kutta code ? Nobody in this forum will be able to understand what we are talking about.

Here it is again for reference:

global dt t_final;
C_p=3.026700000000000e-10;
kp=392400000;
theta_p = 0.231516000000000;
R_s=20*1.8e6; % ohm
ks=1.6e3;
k=(kp*ks)/(kp+ks); 
m=0.1;
xi1=0.000001;
c=2*xi1*sqrt(k*m); 
v_band=2;
uD=0.2; %
uS=0.5; % 
C=20000; % 
Fn=100; % 
dt=1e-5; % time step
t_final=3; % time final
Total_element = ceil(t_final/dt);
% initial condition
t(1)=0;
y1(1)=0;
yy1(1)=v_band;
Vr(1)=0;
frictionf_ary(1)=0;
VCs(1)=0;
for i=1:Total_element
    % update time
    t(i+1)=t(i)+dt;
    
    if rem(i, round(t_final/dt)*0.1)==0    
        disp([' time:', num2str(t(i+1))]);
    end
    
    % calculating vrel
    vrel=v_band-yy1(i); % rel velocity
    if abs(vrel)/v_band < 0.001  % 
        fstatic = uS*Fn*sign(vrel); faply = k*y1(i)+c*yy1(i)+theta_p*(ks/(kp+ks))*Vr(i);
        if abs(faply) < uS*Fn
            frictionf_ary(i+1) = k*y1(i)+c*yy1(i)+theta_p*(ks/(kp+ks))*Vr(i);      
        elseif abs(faply) >= uS*Fn
            frictionf_ary(i+1) = fstatic;             
        end        
    else
        mu=uD+(uS-uD)*exp(-(C*abs(vrel))); % 
        frictionf_ary(i+1)=mu*Fn*sign(vrel);
    end
    
    t_v=t(i); y1_v=y1(i); yy1_v=yy1(i); Vr_v=Vr(i); frictionf=frictionf_ary(i+1); %  
    
    % define function handle
    fy1 = @(t,y1,yy1,Vr) yy1;
    fyy1 = @(t,y1,yy1,Vr) (-c*yy1-k*y1+frictionf-(ks/(kp+ks))*theta_p*Vr)/m;
    fVr = @(t,y1,yy1,Vr) (theta_p*(ks/(kp+ks))*yy1*R_s-Vr)/R_s/C_p;
    
    [k1_fy1, k1_fyy1, k1_fVr, k2_fy1, k2_fyy1, k2_fVr, k3_fy1, k3_fyy1, k3_fVr, k4_fy1, k4_fyy1, k4_fVr] = f_Reg(t_v, y1_v, yy1_v, Vr_v, dt, fy1, fyy1, fVr);
    y1(i+1) = y1(i)+dt/6*(k1_fy1 + 2*k2_fy1 + 2*k3_fy1 + k4_fy1);
    yy1(i+1) = yy1(i)+dt/6*(k1_fyy1 + 2*k2_fyy1 + 2*k3_fyy1 + k4_fyy1);
    Vr(i+1) = Vr(i)+dt/6*(k1_fVr + 2*k2_fVr + 2*k3_fVr + k4_fVr);       
end
%%
%%
figure(2);
plot(t(:),Vr(:),'k'); %voltage
grid on;
xlabel('Time [s]', 'fontsize', 12, 'FontWeight', 'bold', 'fontname', 'Times New Roman');
ylabel('Voltage [V]', 'fontsize', 12, 'FontWeight', 'bold', 'fontname', 'Times New Roman');
%%
%%
end
function [k1_fy1, k1_fyy1, k1_fVr, k2_fy1, k2_fyy1, k2_fVr, k3_fy1, k3_fyy1, k3_fVr, k4_fy1, k4_fyy1, k4_fVr] = f_Reg(t, y1, yy1, Vr, dt, fy1, fyy1, fVr)
k1_fy1  =  fy1(t,           y1,               yy1,                 Vr);
k1_fyy1 =  fyy1(t,          y1,               yy1,                 Vr);
k1_fVr  =  fVr(t,           y1,               yy1,                 Vr);
k2_fy1  =  fy1(t+dt/2,      y1+dt/2*k1_fy1,   yy1+dt/2*k1_fyy1,    Vr+dt/2*k1_fVr);
k2_fyy1  = fyy1(t+dt/2,     y1+dt/2*k1_fy1,   yy1+dt/2*k1_fyy1,    Vr+dt/2*k1_fVr);
k2_fVr  =  fVr(t+dt/2,      y1+dt/2*k1_fy1,   yy1+dt/2*k1_fyy1,    Vr+dt/2*k1_fVr);
k3_fy1  =  fy1(t+dt/2,      y1+dt/2*k2_fy1,   yy1+dt/2*k2_fyy1,    Vr+dt/2*k2_fVr);
k3_fyy1  = fyy1(t+dt/2,     y1+dt/2*k2_fy1,   yy1+dt/2*k2_fyy1,    Vr+dt/2*k2_fVr);
k3_fVr  =  fVr(t+dt/2,      y1+dt/2*k2_fy1,   yy1+dt/2*k2_fyy1,    Vr+dt/2*k2_fVr);
k4_fy1  =  fy1(t+dt,        y1+dt  *k3_fy1,   yy1+dt  *k3_fyy1,    Vr+dt  *k3_fVr);
k4_fyy1  = fyy1(t+dt,       y1+dt  *k3_fy1,   yy1+dt  *k3_fyy1,    Vr+dt  *k3_fVr);
k4_fVr  =  fVr(t+dt,        y1+dt  *k3_fy1,   yy1+dt  *k3_fyy1,    Vr+dt  *k3_fVr);

end
%%

[haohaoxuexi1]

@Torsten

"One error in your code is the line"

"frictionf" depends on all three unknown functions, and you had to include this expression here instead of the constant value "frictionf"."

"It will be quite hard to repair the first one within your Runge-Kutta code form."

I considered your suggestion of changing the expression of frictionf into an expression with "y1 yy1 Vr" instead of a constant, and then it works.

Thanks,

[Torsten]

"frictionf" depends on the result of an if-condition. I can't think of a way to set it correctly in the function handle depending on the result of this if-statement.

[haohaoxuexi1]

@Torsten

"By the way: why did you delete your Runge-Kutta code ? Nobody in this forum will be able to understand what we are talking about."

The code after your modification in you comment included everything to understand the problem. I just don't want to mess thing up.

""frictionf" depends on the result of an if-condition. I can't think of a way to set it correctly in the function handle depending on the result of this if-statement."

In my code, did you notice actually the handle is updated in every iteration, so I just simplily included the handle in the if statement to do it. But yeah, like you said, my code is insufficient and hard to debug.

"I gave you the correct and efficient form of a Runge-Kutta-code of order 4, and the results for your problem were the same as with ODE45."

Your code is more efficient.

[Torsten]

In my code, did you notice actually the handle is updated in every iteration, so I just simplily included the handle in the if statement to do it.

In function "f_Reg", between the calculations of the different ki's, "frictionf" and thus the function handle fyy1 has to be adapted depending on the values of the 2nd, 3rd and 4th argument of the function handles. I don't believe that "simply including the handle in the if statement" can solve this problem.

Another problem is what Jan alluded to: Your if-statements have the effect that your differential equations have discontinuities. Strictly speaking, you would have to stop and restart your integrator in each discontinuity to get reliable results.

[Sam Chak]

+1 Kudos to @Torsten. 👍

[haohaoxuexi1]

@Torsten

Thank you for the useful comment.

For now, my problem is solved because the output is similar to the one using ode45.

I am not an math expert, my purpose is to utilize the tool to solve some problem.

Thank you for giving me this useful advice "frictionf" depends on all three unknown functions, and you had to include this expression here instead of the constant value "frictionf"."

For calculation efficiency related issue, I will try to improve it.

[Sam Chak]

Hi @haohaoxuexi1, getting the RK4 algorithm to yield the same result as the built-in ode45 solver is probably only the first step. Your next step is to ensure that the piecewise-smooth dynamical ODEs are correctly "translated" into the MATLAB code.

I'm also not an expert in many branches of math. However, you have to be clear that good at Math and good at MATLAB are two different things. Good at math, nevertheless, gives you an edge on formulating a mathematical problem in a way that can be efficiently solved by MATLAB.

For this, you are advised to "display" the piecewise-smooth dynamical ODEs in LaTeX so that we can get a clear picture of what exactly the system is. For long equations like yours, the experience is very different from viewing the lines of cluttered codes. Human minds tend to interpret the "display-style" of the Mathematical Equations better than in codes.

[haohaoxuexi1]

@Sam Chak

Thank you.

Answer 2

[Jan]

The function to be integrated contains:

if abs(vrel)/v_band<0.001
    ...
    if abs(faply)<uS*Fn        

This looks like this function is not smooth. Matlab's ODE integrators are designed to handle smooth functions only. Otherwise the step size controller drives mad and the final value is not a "result" anymore but can be dominated by accumulated rounding errors.

A fixed step solver runs brutally over the discontuinuities. It depends on the function and stepsize, if the calculated trajectory is "better" or "worse" than the trajectory calculated with a step size controller. From a scientific point of view, both methods, fixed step RK and ODE45 without detection of discontinuities are expensive random number generators with a weak entropy.

Don't do this for scientific work. Use a stepsize controlled integrator and stop/restart the integration at each discontinuitiy to recondition the controller.

[haohaoxuexi1]

@Jan

"Don't do this for scientific work. Use a stepsize controlled integrator and stop/restart the integration at each discontinuitiy to recondition the controller."

do you have any sample code for this suggestion?

[haohaoxuexi1]

@Jan @Torsten

Hi Jan and Torsten,

I just thought about a new question regarding using ODE45 and rungekutta for solving ODE question.

If I changed Fn=100 to a random number, within 50-150,

F0_0 =100;

F0 = F0_0 + -0.5*F0_0+rand(1,1)*(0.5*F0_0--0.5*F0_0);

Is the method still able to solve the question properly? Or both ode45 and rungekutta won't work well under this situation,

Thanks,

[Torsten]

All deterministic ODE integrators assume that the function to be integrated is at least differentiable. They can't cope with stochastic inputs.

Look up "Stochastic Differential Equations" if you want to define some inputs as random variables.

[haohaoxuexi1]

@Torsten

Thank you for explanation.