# Runga Kutta numerical integration

6 views (last 30 days)
numnum on 20 Oct 2017
Edited: numnum on 20 Oct 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 Comment

the cyclist on 20 Oct 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.

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

### Community Treasure Hunt

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

Start Hunting!

Translated by