ode45 and rungekutta yield different result

1 ビュー (過去 30 日間)
haohaoxuexi1
haohaoxuexi1 2022 年 7 月 10 日
コメント済み: haohaoxuexi1 2022 年 7 月 30 日
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 件のコメント
haohaoxuexi1
haohaoxuexi1 2022 年 7 月 10 日
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
Torsten 2022 年 7 月 10 日
編集済み: Torsten 2022 年 7 月 10 日
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 件のコメント
haohaoxuexi1
haohaoxuexi1 2022 年 7 月 11 日
Thank you.

サインインしてコメントする。

その他の回答 (1 件)

Jan
Jan 2022 年 7 月 10 日
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 件のコメント
haohaoxuexi1
haohaoxuexi1 2022 年 7 月 30 日
Thank you for explanation.

サインインしてコメントする。

Community Treasure Hunt

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

Start Hunting!

Translated by