CODEWITHSUNDEEP
matlabodedifferential-equations
harman singh
harman singh

Reputation: 33

solving system of ode using matlab

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.

Upvotes: 2

Views: 428

Answers (2)

remus
remus

Reputation: 2710

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

Upvotes: 1

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

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

Upvotes: 1

Related Questions