Solving differential equation with Runge Kutta 4th order
15 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Edoardo Bertolotti
le 21 Mar 2022
Commenté : Ahmed J. Abougarair
le 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 commentaires
Réponse acceptée
VBBV
le 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 commentaires
Torsten
le 8 Août 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
le 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)
Plus de réponses (2)
Torsten
le 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 commentaires
Edoardo Bertolotti
le 22 Mar 2022
2 commentaires
Torsten
le 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.
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!