Solving differential equation with Runge Kutta 4th order

32 次查看（过去 30 天）
Edoardo Bertolotti
Edoardo Bertolotti 2022-3-21
评论： Ahmed J. Abougarair 2024-3-24
I've tried to write this code but it seems to give different values than ODE 45.
The equation is the following:
Does it make sense?
Thanks in advance
%Extremes
a=0;
b=1;
%Stepsize
h=0.05;
%time interval and initial value
t=(a:h:b)';
x(1)=0.032;
parameter=0.15;
t1=zeros(size(t));
n=numel(t1);
%RK cycle
for i=1:n-1
%function
dxdt = @(x) parameter*(1-x)-(1-x)*(35*x(i)/((x(i)^2)/0.01+x(i)+1));
j1 = h*dxdt(x(i));
j2 = h*dxdt(x(i)+j1/2);
j3 = h*dxdt(x(i)+j2/2);
j4 = h*dxdt(x(i)+j3);
x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4);
end

请先登录，再回答此问题。

采纳的回答

VBBV
VBBV 2022-3-21
a=0;
b=1;
%Stepsize
h=0.05;
%time interval and initial value
t=(a:h:b)';
x(1)=0.032;
parameter=0.15;
t1=zeros(size(t));
n=numel(t1);
%RK cycle
dxdt = @(x) parameter*(1-x)-(1-x).*(35*x./((x.^2)/0.01+x+1));
for i=1:n-1
%function
j1 = h*dxdt(x(i));
j2 = h*dxdt(x(i)+j1/2);
j3 = h*dxdt(x(i)+j2/2);
j4 = h*dxdt(x(i)+j3);
x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4);
end
plot(t,x)
  4 个评论
显示 2更早的评论
Torsten
Torsten 2022-8-8
Integrate to the first jump point,
set the new initial values to the solution at the last time instant + the jump in the solution,
restart and integrate to the next jump point,
set the new initial values to the solution at the last time instant + the jump in the solution,
...
Ahmed J. Abougarair
Ahmed J. Abougarair 2024-3-24
clc;
clear all;
F = @(t,y) 4*exp(0.8*t)-0.5*y
t0=input('Enter the value of t0 : \n');
y0=input('Enter the value of y0 : \n');
tn=input('Enter the value of t for which you want to find the value of y : \n');
h=input('Enter the step length : \n');
i=0;
while i<tn
k_1 = F(t0,y0);
k_2 = F(t0+0.5*h,y0+0.5*h*k_1);
k_3 = F((t0+0.5*h),(y0+0.5*h*k_2));
k_4 = F(((t0)+h),(y0+k_3*h));
nexty = y0 + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
y0=nexty;
t0=t0+h;
i=i+h;
end
fprintf('The value of y at t=%f is %f \n',t0,y0)
% validate using a decent ODE integrator
tspan = [0,1]; Y0 = 2;
[tx,yx] = ode45(F, tspan, Y0);
fprintf('The true value of y at t=%f is %f \n',tspan(end),yx(end))
Et= (abs(yx(end)-y0)/yx(end))*100;
fprintf('The value of error Et at t=%f is %f%% \n',tspan(end),Et)

请先登录，再进行评论。

更多回答（2 个）

Torsten
Torsten 2022-3-21
y0 = 0.3075;
tspan = [0 6];
[t,y] = ode45(@trial_ODE45,tspan,y0);
plot(t,y,'-')
function dydt = trial_ODE45(t,y)
parameter = 1.5;
%equations system
dydt = parameter*(1-y)-(1-y).*(35*y./((y.^2)/0.01+y+1));
end
Edoardo Bertolotti
Edoardo Bertolotti 2022-3-22
Ok, thanks to both, the code is working now, but still, like the main problem, it shows the same results:
Euler is similar to ODE45, but RK4 no..I would expect ODE45 to be more similar to RK than Euler or equal to both, using same steps etc
a=0;
b=6;
%Stepsize
h=0.05;
%time interval and initial value
t=(a:h:b)';
x(1)=0.3075;
parameter=1.5;
t1=zeros(size(t));
n=numel(t1);
%Euler
dxdt = @(x) parameter*(1-x)-(1-x).*(35*x./((x.^2)/0.01+x+1));
for i=1:n-1
x(i+1)=x(i)+h*dxdt(x(i));
end
ylim([-0.1 1.1])
hold on
plot(t,x)
title('methods')
hold on
grid on
xlabel('time')
ylabel('concentration')
hold on
%RK cycle
dxdt = @(x) parameter*(1-x)-(1-x).*(35*x./((x.^2)/0.01+x+1));
for i=1:n-1
%function
j1 = h*dxdt(x(i));
j2 = h*dxdt(x(i)+j1/2);
j3 = h*dxdt(x(i)+j2/2);
j4 = h*dxdt(x(i)+j3);
x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4);
end
hold on
plot(t,x)
hold on
%ODE45
y0 = 0.3075;
tspan = [0 6];
[t,y] = ode45(@trial_ODE45,tspan,y0);
plot(t,y,'-')
legend('Euler','Runge Kutta','ODE45','Location','NorthOutside','Orientation','horizontal','Box','off')
% % %runge kutta
function dydt = trial_ODE45(t,y)
parameter = 1.5;
%equations system
dydt = parameter*(1-y)-(1-y).*(35*y./((y.^2)/0.01+y+1));
end

类别

在 Help Center 和 File Exchange 中查找有关 Runge Kutta Methods 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by