# ode45 and rungekutta yield different result

1 view (last 30 days)
haohaoxuexi1 on 10 Jul 2022
Commented: haohaoxuexi1 on 30 Jul 2022
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.
Thanks,
haohaoxuexi1 on 10 Jul 2022
I didn't get what you mean, can you please show me a simple example how to do this?
I am new to this.

Torsten on 10 Jul 2022
Edited: Torsten on 10 Jul 2022
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
haohaoxuexi1 on 11 Jul 2022
Thank you.

Jan on 10 Jul 2022
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 on 30 Jul 2022
Thank you for explanation.

### Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

R2018b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!