# ode45 and rungekutta yield different result

1 ビュー (過去 30 日間)
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.
Thanks,
##### 5 件のコメント表示非表示 4 件の古いコメント
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 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 件のコメント表示非表示 11 件の古いコメント
haohaoxuexi1 2022 年 7 月 11 日
Thank you.

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

### その他の回答 (1 件)

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 件のコメント表示非表示 3 件の古いコメント
haohaoxuexi1 2022 年 7 月 30 日
Thank you for explanation.

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

### カテゴリ

Find more on Structural Mechanics 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!

Translated by