solving system of ode using matlab -
i have 9 equations time dependent coefficient g
% 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 , g = 0.1*exp(-((t-200)/90).^2) time dependent coefficient , time t = 0:500; initial condition [0 0 1 0 0 0 0 0 0].
i'm getting wrong negative values output y(1), y(2). can pls try solve above eqns ode45 can compare results.
with simple application of rk4 picture
nicely positive, 1 strange initial jump in y(1) component. note scale, on whole y(1) rather small. seems system stiff @ point, rk45 might have problems, implicit runge-kutta method better.
and zoom of initial oscillations
python code
import numpy np import matplotlib.pyplot plt 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] t0 in t[0:-1]: n=200; h = dt/n; 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() 
Comments
Post a Comment