runge-kutta ode45

37 vues (au cours des 30 derniers jours)
Alvaro Mª Zumalacarregui Delgado
I want to solve a two EDO system with the runge kutta method, I wrote this code but i don't get anything, someone know why this is wrong? if i use command ode45 is better
function [] = runge_kutta_sis(f,g,x0,y0,a,b,h)
t = a:h:b; n = length (t);
x = (x0); y = (y0);
for i=1:n-1
k1 = h*f(x(i),y(i),t(i));
l1 = h*g(x(i),y(i),t(i));
k2 = h*f(x(i)+k1/2,y(i)+l1/2,t(i)+h/2);
l2 = h*g(x(i)+k1/2,y(i)+l1/2,t(i)+h/2);
k3 = h*f(x(i)+k2/2,y(i)+l2/2,t(i)+h/2);
l3 = h*g(x(i)+k2/2,y(i)+l2/2,t(i)+h/2);
k4 = h*f(x(i)+k3,y(i)+l3,t(i)+h);
l4 = h*g(x(i)+k3,y(i)+l3,t(i)+h);
x (i+1) = x(i) + (1/6)*(k1+2*k2+2*k3+k4);
y (i+1) = y(i) + (1/6)*(l1+2*l2+2*l3+l4);
end
plot (t,x,t,y)
I write this in the command window but I can't see any figure
>> f = @(x,y,t) 10-0.1*x*y;
>> g = @(x,y,t) 20-0.2*x*y;
>> runge_kutta_sis(f,g,500,600,0,20,0.5)

Réponse acceptée

KALYAN ACHARJYA
KALYAN ACHARJYA le 1 Mar 2021
Minor modification: Function Code (Save in the same working directory, filename runge_kutta_sis.m)
function runge_kutta_sis(f,g,x0,y0,a,b,h)
t=a:h:b;
n=length(t);
x=zeros(1,length(n));
y=zeros(1,length(n));
x(1)=x0;
y(1)=y0;
for i=1:n-1
k1 = h*f(x(i),y(i),t(i));
l1 = h*g(x(i),y(i),t(i));
k2 = h*f(x(i)+k1/2,y(i)+l1/2,t(i)+h/2);
l2 = h*g(x(i)+k1/2,y(i)+l1/2,t(i)+h/2);
k3 = h*f(x(i)+k2/2,y(i)+l2/2,t(i)+h/2);
l3 = h*g(x(i)+k2/2,y(i)+l2/2,t(i)+h/2);
k4 = h*f(x(i)+k3,y(i)+l3,t(i)+h);
l4 = h*g(x(i)+k3,y(i)+l3,t(i)+h);
x(i+1)=x(i)+(1/6)*(k1+2*k2+2*k3+k4);
y(i+1)=y(i)+(1/6)*(l1+2*l2+2*l3+l4);
end
plot(t,x,t,y);
grid on;
Main Script:
f = @(x,y,t) 10-0.1*x*y;
g = @(x,y,t) 20-0.2*x*y;
runge_kutta_sis(f,g,500,600,0,20,0.5)
Output:
Hope it Helps!
  2 commentaires
Alvaro Mª Zumalacarregui Delgado
it helps thanks! if I want to do it with the command ode45 how can I write it?
KALYAN ACHARJYA
KALYAN ACHARJYA le 1 Mar 2021

Connectez-vous pour commenter.

Plus de réponses (1)

Meysam Mahooti
Meysam Mahooti le 5 Mai 2021
https://www.mathworks.com/matlabcentral/fileexchange/61130-runge-kutta-fehlberg-rkf78?s_tid=srchtitle

Community Treasure Hunt

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

Start Hunting!

Translated by