ode45 - Di erential Equation Solver - Purdue University

1 ode45 - di erential Equation SolverThis routine uses a variable step Runge-Kutta Method to solve di erential equations syntax forode45for rst order di erential equations and that for second order di erentialequations are basically the same. However, the .m les are quite di Order equations (y0=f(t;y)y(t0)=y0A. Create a .m le forf(t;y) (see the tutorial on numerical methods and m les on how todo this). Save le as, for example, Basic syntax type :[t,y]= ode45 ('yp',[t0,tf],y0);(your version ofode45may not require brackets around t0, tf)8> <>:yp = the .m file of the functionf(t;y)saved as , tf = initial and terminal values ofty0 = initial value ofyatt0C. For example, to numerically solve(t2y0=y+3ty(1) = 2over the interval 1 t 4: Create and save the le for the function1t2(y+3t). At aMatlabprompt type:[t,y]= ode45 ('yp',[1,4],-2);(your version ofode45may not require backets around [1,4]) To print results type :[t,y] To plot results type :plot(t,y) To plot results type with a '+' symbol :plot(t,y,'+') Order Equations8> <>:y00+p(t)y0+q(t)y=g(t)y(t0)=y0y0(t0)=y1A .))

2 First convert 2ndorder Equation above to an equivalent system of 1storder ; x2=y0:(x01=x2x02= q(t)x1 p(t)x2+g(t);wherex1(t0)=y0;x2(t0)=y1:B. Create and save a .m le which will return a vector-valuedfunction. This is a little trickyso here is a speci c example. Suppose the system is as below and 0 t 4(x01=x2x02= tx1+etx2+3sin2t;wherex1(0)=2;x2(0) = 8: Create the following function le and save it :function xp=F(t,x)xp=zeros(2,1); % since output must be a column vectorxp(1)=x(2);xp(2)=-t*x(1)+exp(t)*x( 2)+3*sin(2*t); % don't forget ; after each line Basic syntax , type :[t,x]= ode45 ('F',[t0,tf],[x10,x20]);8>>> <>>>:F = the .m file of the vector-function saved as abovet0, tf = initial and terminal values oftx10 = initial value ofx1att0: x10 =x1(t0)x20 = initial value ofx2att0: x20 =x2(t0)The example above becomes :[t,x]= ode45 ('F',[0,4],[2,8]); Sincex1(t)=y, to print out the values of the solutionyfort0 t tf,type:[t,x(:,1)]To plot the solution on a graphtvsy,type:plot(t,x(:,1))(This is because the vectorxhas 1stcomponentx1=yand 2ndcomponentx2=y0.)))

3 To p l o tx1vsx2(phase plane) type :plot(x(:,1),x(:,2))