Solving differential equation with Runge Kutta 4th order
24 views (last 30 days)
Show older comments
Edoardo Bertolotti
on 21 Mar 2022
Commented: Ahmed J. Abougarair
on 24 Mar 2024
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
0 Comments
Accepted Answer
VBBV
on 21 Mar 2022
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 Comments
Torsten
on 8 Aug 2022
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
on 24 Mar 2024
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)
More Answers (2)
Torsten
on 21 Mar 2022
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
0 Comments
Edoardo Bertolotti
on 22 Mar 2022
2 Comments
Torsten
on 22 Mar 2022
x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3+j4);
instead of
x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4);
in RK4.
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!