Runga Kutta numerical integration

4 visualizzazioni (ultimi 30 giorni)
numnum
numnum il 20 Ott 2017
Modificato: numnum il 20 Ott 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 Commento
the cyclist
the cyclist il 20 Ott 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.

Accedi per commentare.

Risposta accettata

Mischa Kim
Mischa Kim il 20 Ott 2017
numnum, check out this answer for the Lorentz attractor.

Più risposte (0)

Community Treasure Hunt

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

Start Hunting!

Translated by