602 CHAPTER 31. NUMERICAL SOLUTIONS FOR SYSTEMS

using the improved Euler method, do the following:

k1 ≡ f (tk,yk) ,k2 ≡ f (tk +h,yk +k1h) , tk = kh, yk+1 = yk +h2(k1+k2)

then there is some constant C such that |y (kh)−yk|<Ch2.

It predicts what the slope should be at a point and then averages the two values toadvance another step.

This problem of getting solutions to first order systems of differential equations hasbeen studied extensively and a book like this is not the place to see a careful description ofthe best methods. However, one of the very best was developed long before computers byRunge and Kutta in 1901.

PROCEDURE 31.1.2 To find a numerical solution to

y′ = f (t,y) , y (0) = y0

using the Runge-Kutta algorithm, do the following: For tk = kh,k = 0, · · · ,

k1 ≡ f (tk,yk) ,k2 ≡ f

(tk +

h2,yk +k1

h2

),

k3 ≡ f

(tk +

h2,yk +k2

h2

),k4 ≡ f (tk +h,yk +k3h)

yk+1 = yk +h6(k1 +2k2 +2k3 +k4)

then there is some constant C such that |y (kh)−yk|<Ch4.

You can have MATLAB use the Runge-Kutta algorithm to numerically find a solution toa system of ordinary differential equations. I will illustrate with a first order system whichcomes from the Van der Pol equation. The exact equation studied is not too important atthis point. My intent is to illustrate the syntax used. Here it is:

f=@(t,x)[x(2),-((x(1)ˆ2-1)*x(2)+x(1))]; n=300; h=.05;y(1,:)=[1,0]; t(1)=0;hold onfor r=1:nk1=f(t(r),y(r,:)); k2=f(t(r)+h/2,y(r,:)+k1*(h/2));k3=f(t(r)+h/2,y(r,:)+k2*(h/2)); k4=f(t(r)+h,y(r,:)+k3*h);y(r+1,:)=y(r,:)+(h/6)*(k1+2*k2+2*k3+k4); t(r+1)=t(r)+h;endplot(t,y(:,1),t,y(:,2))interp1(t,[y(:,1),y(:,2)],[2,3,4,5,6])

Then press “enter”. This will compute the solution to the differential equation(xy

)′=

(y

−((

x2−1)

y+ x) ) ,

(x(0)y(0)

)=

(10

).