Runga Kutta numerical integration

5 ビュー (過去 30 日間)
numnum
numnum 2017 年 10 月 20 日
編集済み: numnum 2017 年 10 月 20 日
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 件)

カテゴリ

Help Center および File ExchangeNumerical Integration and Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by