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
– Lutz Lehmann, Commented Apr 1, 2015 at 8:07
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
– harman singh, Commented Apr 1, 2015 at 8:55
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
– Lutz Lehmann, Commented Apr 1, 2015 at 9:13
dear Lutzl , Thanks ! your plots seems in agreement with the results i am looking for .do you used Odeint in python ? — harman singh
– harman singh, Commented Apr 1, 2015 at 10:58
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
– remus, Commented Apr 1, 2015 at 11:45

Lutz Lehmann · Accepted Answer · 2015-04-01 11:24:22Z

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()

remus · Accepted Answer · 2015-04-01 11:43:43Z

1

And in Matlab:

options = odeset('AbsTol', 1e-12);
[T,Y] = ode45(@tarak, [0, 500], [0 0 1 0 0 0 0 0 0], options);

enter image description here

answered Apr 1, 2015 at 11:43

remus

2,7102 gold badges25 silver badges52 bronze badges

Collectives™ on Stack Overflow

solving system of ode using matlab

2 Answers 2

Python code

Comments

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Python code

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related