ode45 and rungekutta yield different result

6 views (last 30 days)
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,
  5 Comments
Torsten
Torsten on 10 Jul 2022
Edited: Torsten on 10 Jul 2022
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
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.

Sign in to comment.

Accepted Answer

Torsten
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
  12 Comments
Sam Chak
Sam Chak on 11 Jul 2022
Edited: Sam Chak on 11 Jul 2022
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.

Sign in to comment.

More Answers (1)

Jan
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.
  4 Comments
Torsten
Torsten on 29 Jul 2022
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.

Sign in to comment.

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!