Runga Kutta numerical integration

5 visualizaciones (últimos 30 días)
numnum
numnum el 20 de Oct. de 2017
Editada: numnum el 20 de Oct. de 2017
Hi
I am trying to use the Runga Kutta method to solve 3 differential equations:
_dg/dt = [-g + f(i_c - i -h)]/b, where f is a function f(x)=0 if x<0, f(x)=x if 0<=x<a, f(x)=a if vm<=x
_dh/dt= (f*g -h)/d
_di/dt=(e*g -i)/c
where i_c, b, d,f,e and c are constants.
But this not yielding the expected result (a rapid rise in g when i_c changes from 1.15 to 3.1 followed by a rapid exponential decrease, and then a slow exponentatial decrease is expected) Does anyone know what might be wrong?
%clear
clc;
clear;
%constants
a = 4;
b = 0.2;
c = 24;
d=294;
e=6;
f=9.4;
k=1.15;
l=3.1;
%function handles
F_g=@(g,i_c,h,i,t) ((i_c -h -i)<0)*(-g/b)+ ((i_c -h -i)>=0)*((i_c -h -i)<a)*((i_c -h -i - g)/b) + ((i_c -h -i)>=a)*((a -g)/b);
F_h=@(h,g,t) (f*g -h)/d;
F_i=@(i,g,t) (e*g -i)/c;
%step
h=0.01;
t = 0:h:1000;
%prealocate memory
g=zeros(1,length(t));
h=zeros(1,length(t));
i=zeros(1,length(t));
i_c=zeros(1,length(t));
k_1_g=zeros(1,length(t));
k_2_g=zeros(3,length(t));
k_3_g=zeros(3,length(t));
k_4_g=zeros(3,length(t));
k_1_h=zeros(1,length(t));
k_2_h=zeros(2,length(t));
k_3_h=zeros(2,length(t));
k_4_h=zeros(2,length(t));
k_1_i=zeros(1,length(t));
k_2_i=zeros(2,length(t));
k_3_i=zeros(2,length(t));
k_4_i=zeros(2,length(t));
%initial values
g(1)=k/(1+e+f);
i(1)=e*g(1);
h(1)=f*g(1);
i_c(1:3000)=l;
i_c(3001:5000)=k;
i_c(5001:length(t))=l;
for idx=2:(length(t))
k_1_g(idx) = F_g(t(idx-1),g(idx-1),h(idx-1),idx(idx-1));%v
k_1_h(idx) = F_h(t(idx-1),g(idx-1),h(idx-1));%as
k_1_i(idx) = F_i(t(idx-1),g(idx-1),i(idx-1));%af
k_2_g(idx) = F_g(t(idx-1)+0.5*h,g(idx-1)+0.5*h*k_1_g(idx),h(idx-1)+0.5*h*k_1_h(idx),i(idx-1)+0.5*h*k_1_i(idx));
k_2_h(idx) = F_h(t(idx-1)+0.5*h,h(idx-1)+0.5*h*k_1_h(idx),g(idx-1)+0.5*h*k_1_g(idx));
k_2_i(idx) = F_i(t(idx-1)+0.5*h,idx(idx-1)+0.5*h*k_1_i(idx),g(idx-1)+0.5*h*k_1_g(idx));
k_3_g(idx) = F_g((t(idx-1)+0.5*h),(g(idx-1)+0.5*h*k_2_g(idx)),(h(idx-1)+0.5*h*k_2_h(idx)),(i(idx-1)+0.5*h*k_2_i(idx)));
k_3_h(idx) = F_h((t(idx-1)+0.5*h),(h(idx-1)+0.5*h*k_2_h(idx)),(g(idx-1)+0.5*h*k_2_g(idx)));
k_3_i(idx) = F_i((t(idx-1)+0.5*h),(i(idx)+0.5*h*k_2_h(idx)),(g(idx-1)+0.5*h*k_2_g(idx)));
k_4_g(idx) = F_g((t(idx-1)+h),(g(idx-1)+k_3_g(idx)*h), (h(idx-1)+k_3_h(idx)*h), (i(idx-1)+k_3_i(idx)*h));
k_4_h(idx) = F_h((t(idx-1)+h),(h(idx-1)+k_3_h(idx)*h), (g(idx-1)+k_3_g(idx)*h));
k_4_i(idx) = F_i((t(idx-1)+h),(i(idx-1)+k_3_h(idx)*h),(g(idx-1)+k_3_g(idx)*h));
g(idx) = g(idx-1) + (1/6)*(k_1_g(idx)+2*k_2_g(idx)+2*k_3_g(idx)+k_4_g(idx))*h;
h(idx) = h(idx-1) + (1/6)*(k_1_h(idx)+2*k_2_h(idx)+2*k_3_h(idx)+k_4_h(idx))*h;
i(idx) = i(idx-1) + (1/6)*(k_1_i(idx)+2*k_2_i(idx)+2*k_3_i(idx)+k_4_i(idx))*h;% main equation
end
plot(t,g)
  1 comentario
the cyclist
the cyclist el 20 de Oct. de 2017
Your code won't execute for me, because of the problem that this function
F_i=@(i,v,t) (e*g -i)/c;
does not have a sensible definition. I think you've messed up your input definitions, so it is looking for a constant g.

Iniciar sesión para comentar.

Respuesta aceptada

Mischa Kim
Mischa Kim el 20 de Oct. de 2017
numnum, check out this answer for the Lorentz attractor.

Más respuestas (0)

Categorías

Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by