Solving differential equation with Runge Kutta 4th order
15 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Edoardo Bertolotti
el 21 de Mzo. de 2022
Comentada: Ahmed J. Abougarair
el 24 de Mzo. de 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 comentarios
Respuesta aceptada
VBBV
el 21 de Mzo. de 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 comentarios
Torsten
el 8 de Ag. de 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
el 24 de Mzo. de 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)
Más respuestas (2)
Torsten
el 21 de Mzo. de 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 comentarios
Edoardo Bertolotti
el 22 de Mzo. de 2022
2 comentarios
Torsten
el 22 de Mzo. de 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.
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!