11—Numerical Analysis 27711.5 Differential Equations
To solve the first order differential equation
y′=f(x,y) y(x 0 ) =y 0 , (11.31)
the simplest algorithm is Euler’s method. The initial conditions arey(x 0 ) =y 0 , andy′(x 0 ) =f(x 0 ,y 0 ),
and a straight line extrapolation is
y(x 0 +h) =y 0 +hf(x 0 ,y 0 ). (11.32)
You can now iterate this procedure using this newly found value ofyas a new starting condition to go
fromx 0 +htox 0 + 2h.
Runge-Kutta
Euler’s method is not very accurate. For an improvement, change from a straight line extrapolation
to a parabolic one. Takex 0 = 0to keep the algebra down, and try a solution near 0 in the form
y(x) =α+βx+γx^2 ; evaluateα,β, andγso that the differential equation is satisfied nearx= 0,
y′=β+ 2γx=f(x,α+βx+γx^2 ).
Recall the Taylor series expansion for a function of two variables, section2.5:
f(x,y) =f(x 0 ,y 0 ) + (x−x 0 )D 1 f(x 0 ,y 0 )+(y−y 0 )D 2 f(x 0 ,y 0 ) +
1
2
(x−x 0 )^2 D 1 D 1 f(x 0 ,y 0 )
+
1
2
(y−y 0 )^2 D 2 D 2 f(x 0 ,y 0 ) + (x−x 0 )(y−y 0 )D 1 D 2 f(x 0 ,y 0 ) +··· (11.33)
=⇒ β+ 2γx=f(0,α) +xD 1 f(0,α) + (βx+γx^2 )D 2 f(0,α) +···. (11.34)
The initial condition is atx= 0,y=y 0 , soα=y 0. Equate coefficients of powers ofxas high as is
possible (here throughx^1 ).
β=f(0,α) 2 γ=D 1 f(0,α) +βD 2 f(0,α).
(If you setγ= 0, this is Euler’s method.)
y(h) =y 0 +hf(0,y 0 ) +
h^2
2
[
D 1 f(0,y 0 ) +f(0,y 0 )D 2 f(0,y 0 )
]
. (11.35)
The next problem is to evaluate these derivatives. If you can easily do them analytically, you canchoose to do that. Otherwise, since they appear in a term that is multiplied byh^2 , it is enough to use
the simplest approximation for the numerical derivative,
D 1 f(0,y 0 ) =
[
f(h,y 0 )−f(0,y 0 )
]
/h (11.36)
You cannot expect to use the same interval,h, for theyvariable — it might not even have the same
dimensions,