Mathematical Tools for Physics

11—Numerical Analysis 335

y′=yis
y(0) =− 3. 6 y(−h) + 5. 2 y(− 2 h),

or in terms of an index notation
yn=− 3. 6 yn− 1 + 5. 2 yn− 2.

This is a linear, constant coefficient, difference equation, and the method for solving it is essentially the same as
for a linear differential equation — assume an exponential formyn=kn.

kn=− 3. 6 kn−^1 + 5. 2 kn−^2

k^2 + 3. 6 k− 5 .2 = 0

k= 1.11 and − 4. 71

Just as with the differential equation, the general solution is a linear combination of these two functions ofn:

yn=A(1.11)n+B(− 4 .71)n,

whereAandBare determined by two conditions, typically specifyingy 1 andy 2. IfB= 0, thenynis proportional
to 1. 11 nand it is the well behaved exponential solution that you expect. If, however, there is even a little bit ofB
present (perhaps because of roundoff errors), that term will eventually dominate and cause the large oscillations.
IfBis as small as 10 −^6 , then whenn= 9the unwanted term is greater than 1.
When I worked out the coefficients in Eq. ( 42 ) the manipulations didn’t look all that different from those
leading to numerical derivatives or integrals, but the result was useless. This is a caution. You’re in treacherous
territory here; tread cautiously.
Are Adams-type methods useless? No, but you have to modify the development in order to get a stable
algorithm. The difficulty in assuming the form

y(0) =

∑N

1

αky(−kh) +

∑N

1

βky′(−kh)

is that the coefficientsαkare too large. To cure this, you can give up some of the 2 Ndegrees of freedom that
the method started with, and pick theαka priorito avoid instability. There are two common ways to do this,