solving system of ode using matlab

Question

I have 9 equations with a time dependent coefficient g

% MY M file
function dy =tarak(t,y)
G= 3.16;
g =  0.1*exp(-((t-200)/90).^2);
dy=zeros(9,1);
dy(1)=-2*2*y(1)+2*G*y(5)+2*g*y(7);
dy(2)=2*y(1)-2*G*y(5);
dy(3)=2*y(1)-2*g*y(7);
dy(4)=-2*y(4)+g*y(9);
dy(5)=-2*y(5)+G*(y(2)-y(1))+g*y(8);
dy(6)=-2*y(6)-G*y(9);
dy(7)=-2*y(7)+g*(y(3)-y(1))+G*y(8);
dy(8)=-G*y(7)-g*y(5);
dy(9)=G*y(6)-g*y(4);

then in command window:

[T,Y] = ode45(@tarak,[0 ,500],[0 0 1 0 0 0 0 0 0])

where coefficient G = 3.16 and g = 0.1*exp(-((t-200)/90).^2) is a time dependent coefficient and time t = 0:500; Initial condition [0 0 1 0 0 0 0 0 0].

I'm getting WRONG negative values for output y(1), y(2). Can someone pls try to solve above eqns with ode45 so that i can compare the results.

Where do you get that these components have to stay positive? Does that not depend on y(5) and y(7)? — Lutz Lehmann
hi Lutzl , Y(1) ,Y(2), Y(3) are probabilities hence can't be negative , my friend solved above equations in C and getting it right . you too are getting negative value ? — harman singh
No, see answer. The code was in python, but that should not matter. You should post a more complete code to see if there are any easily discernable problems. — Lutz Lehmann
dear Lutzl , Thanks ! your plots seems in agreement with the results i am looking for .do you used Odeint in python ? — harman singh
Matlab code is ok but you need to play with the tolerances because your first state has an order of 1e-8 and the default abs tolerance of ode45 is 1e-6. — remus

Lutz Lehmann Lutz Lehmann · Accepted Answer · 2015-04-01T08:54:22

With a simple application of RK4 I get this picture

enter image description here

nicely positive, with one strange initial jump in the y(1) component. But note the scale, on the whole y(1) is rather small. It seems that the system is stiff at this point, so rk45 might have problems, an implicit Runge-Kutta method would be better.

And a zoom of the initial oscillations

enter image description here

Python code

import numpy as np
import matplotlib.pyplot as plt
from math import exp

def dydt(t,y):
    dy = np.array(y);

    G = 3.16;
    g = 0.1*exp(-((t-200)/90)**2);

    dy[0]=-2*2*y[0]+2*G*y[4]+2*g*y[6];
    dy[1]=   2*y[0]-2*G*y[4];
    dy[2]=   2*y[0]-2*g*y[6];
    dy[3]=  -2*y[3]+  g*y[8];
    dy[4]=  -2*y[4]+  G*(y[1]-y[0])+g*y[7];
    dy[5]=  -2*y[5]-  G*y[8];
    dy[6]=  -2*y[6]+  g*(y[2]-y[0])+G*y[7];
    dy[7]=  -G*y[6]-  g*y[4];
    dy[8]=   G*y[5]-  g*y[3];
    return dy;

def RK4Step(f,x,y,h):
    k1=f(x      , y         )
    k2=f(x+0.5*h, y+0.5*h*k1)
    k3=f(x+0.5*h, y+0.5*h*k2)
    k4=f(x+    h, y+    h*k3)
    return (k1+2*(k2+k3)+k4)/6.0


t= np.linspace(0,500,200+1);
dt = t[1]-t[0];
y0=np.array([0, 0, 1, 0, 0, 0, 0, 0, 0]);

y = [y0]

for t0 in t[0:-1]:
    N=200;
    h = dt/N;
    for i in range(N):
        y0 = y0 + h*RK4Step(dydt,t0+i*h,y0,h);
    y.append(y0);

y = np.array(y);

plt.subplot(121);
plt.title("y(1)")
plt.plot(t,y[:,0],"b.--")
plt.subplot(122);
plt.title("y(2)")
plt.plot(t,y[:,1],"b-..")
plt.show()

solving system of ode using matlab

2 Answers

Python code