Runga Kutta numerical integration

4 次查看(过去 30 天)
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 个评论
the cyclist
the cyclist 2017-10-20
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.

请先登录,再进行评论。

采纳的回答

Mischa Kim
Mischa Kim 2017-10-20
numnum, check out this answer for the Lorentz attractor.

更多回答(0 个)

Community Treasure Hunt

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

Start Hunting!

Translated by