by R. Grothmann

We investigate the stability of the BDF (backward differential formula) following the work of Gear e.a.

A BDF is a formula to solve differential equation. It uses multiple steps, and it is implicit. We assume the solution in the points x_0,...,x_n given as y_0,...,y_n and determine y_{n+1} implicitely such that

where p interpolations (x_i,y_i) for i=0,...n, and satisfies

First, let us determine the formula for p for the case n=1.

>function p(x) &= a+b*x+c*x^2;

It suffices to solve the special case

since we are only interested in equidistant discretization here.

>sol &= solve([p(-h)=y0,p(0)=y1,diffat(p(x),x=h)=yd2],[a,b,c])

h yd2 + 2 y1 - 2 y0 h yd2 - y1 + y0 [[a = y1, b = -------------------, c = ---------------]] 3 h 2 3 h

We define a function for the solution.

>function ps(x) &= p(x) with sol[1]

2 x (h yd2 - y1 + y0) x (- h yd2 - 2 y1 + 2 y0) -------------------- - ------------------------- + y1 2 3 h 3 h

The value of this function at h defines p(x_{n+1}) from the previous values and the derivative p'(x_{n+1}).

>function ystep(y0,y1,yd2,h) &= ps(h)|ratsimp

2 h yd2 + 4 y1 - y0 ------------------- 3

Here is an Euler function, which solves differential equations using this two step BDF. We use iteration for the solution of the non-linear equation

This is not optimal and must be replaced with faster solvers in real life.

>function bdf (f,a,b,n,y0) ... x=linspace(a,b,n); y=zeros(size(x)); y[1:2]=runge(f,x[1:2],y0); h=(b-a)/n; for i=2 to n; y[i+1]=y[i]+f(x[i],y[i])*h; repeat yold=y[i+1]; yd2=f(x[i+1],yold); y[i+1]=&:ystep(y[i-1],y[i],yd2,h); until yold~=y[i+1]; end; end return y; endfunction

Let us show that the order of this procedure is 2.

>error1000=bdf("-2*x*y",0,1,1000,1)[-1]-1/E

4.9057824697e-007

>error2000=bdf("-2*x*y",0,1,2000,1)[-1]-1/E

1.22635597177e-007

>log(error1000/error2000)/log(2)

2.00010545601

Next, we solve the special equation

The non-linear problem for p can be solved explicitely for this equation. We get an explicit equation for y1 given y_0,y_1 to y_2.

>difgl &= solve(ystep(y0,y1,lambda*y2,h)=y2,y2)

y0 - 4 y1 [y2 = --------------] 2 h lambda - 3

Let us show for lambda=-5, that this recursion does indeed converge to the correct solution.

>lambda=-5; n=1000; h=1/n; ... sequence("(x[n-2]-4*x[n-1])/(2*h*lambda-3)",[1,exp(lambda/n)],n+1)[-1]

0.00673766561896

>exp(-5)

0.00673794699909

For stability, a formula must be able to solve this case, even if h*lambda is large. We want the iteration to converge to 0 for each lambda<0, or even better, for each complex lambda with re(lambda)<0. This is called A-stability.

The two step PDF is stable for this lambda and h, as we see with the following start values.

>sequence("(x[n-2]-4*x[n-1])/(2*h*lambda-3)",[1,0],n+1)[-1]

-0.00340277570995

>sequence("(x[n-2]-4*x[n-1])/(2*h*lambda-3)",[0,1],n+1)[-1]

0.0101912705026

We can study the problem exactly. The sequence of y_n is defined by a so-called difference formula. To solve such formulas, we can make the ansatz y_n=t^n. It follows easily, that solutions of this formula are determined by the zeros of

If t is a zero, t^n will satisfy the difference formula. We are interested to see |t|<1, if re(lambda*h)<0.

Replacing t=1/z, we seek the solutions of

We want to show that no |z| <= 1 with re(p(z))<0 exists.

>function q(z) &= 3-4*z+z^2

2 z - 4 z + 3

Let us plot the image of |z|=1. We see that indeed, the method is A-stable, since the interior of the image is in the right half plane.

>z=exp(I*linspace(0,2pi,1000)); plot2d(q(z),r=10):

Discussing that involves a discussion of the function

>&realpart(q(cos(x)+I*sin(x)))|trigsimp

2 2 - sin (x) + cos (x) - 4 cos(x) + 3

We can do that analytically, but here is a simple plot.

>plot2d("2*x^2-4*x+2",-1,1):

Next we do the same for the three term BDF.

>function p(x) &= a+b*x+c*x^2+d*x^3; >sol &= solve([p(-2*h)=y0,p(-h)=y1,p(0)=y2,diffat(p(x),x=h)=yd3],[a,b,c,d]); >function ps(x) &= p(x) with sol[1]; >function ystep(y0,y1,y2,yd3,h) &= ps(h)|ratsimp

6 h yd3 + 18 y2 - 9 y1 + 2 y0 ----------------------------- 11

>sol &= solve(ystep(y0,y1,y2,lambda*y3,h)=y3,y3)

- 18 y2 + 9 y1 - 2 y0 [y3 = ---------------------] 6 h lambda - 11

The formula for

has three terms now. We follow the same ideas as above.

>function q(z) &= 11+(num(rhs(sol[1])) with [y0=z^3,y1=z^2,y2=z])

3 2 - 2 z + 9 z - 18 z + 11

This polynomial now stretches into the left half plane. This makes the method unstable for h*lambda=i.

>plot2d(q(z),r=50):

We derive an instable and unusable formula. The idea is to interpolate

using Hermite interpolation, then to set

>function p(x) &= a+b*x+c*x^2+d*x^3

3 2 d x + c x + b x + a

>sol &= solve( ... [p(-h)=y0,diffat(p(x),x=-h)=yd0,p(0)=y1,diffat(p(x),x=0)=yd1],[a,b,c,d])

h (2 yd1 + yd0) - 3 y1 + 3 y0 [[a = y1, b = yd1, c = -----------------------------, 2 h h (yd1 + yd0) - 2 y1 + 2 y0 d = ---------------------------]] 3 h

We can easily derive a formula for the next y-value.

>function ystep(y0,yd0,y1,yd1,h) &= p(h) with sol[1]

h (2 yd1 + yd0) + h (yd1 + yd0) + h yd1 - 4 y1 + 5 y0

To investigate the stability, we use y'=lambda*y as above.

>y2 &= ystep(y0,lambda*y0,y1,lambda*y1,h)

h (2 y1 lambda + y0 lambda) + h (y1 lambda + y0 lambda) + h y1 lambda - 4 y1 + 5 y0

Set h*lambda=c.

>y2c &= ratsimp(y2 with lambda=c/h)

(4 c - 4) y1 + (2 c + 5) y0

This difference recursion should be stable for re(c)<0. So let us investigate this.

>&c with solve(y2n-y2c,c), function h(t) &= % with [y2n=t^2,y1=t,y0=1]

y2n + 4 y1 - 5 y0 ----------------- 4 y1 + 2 y0 2 t + 4 t - 5 ------------ 4 t + 2

We now do not want to have values |x| > 1, which touch the negative axis.

>z=exp(I*linspace(0,2pi,1000)); ... plot2d(h(z),r=6):

But we see that such x exist for all c with re(c)<0. Indeed, the algorithm fails completely, even for moderate c.

>lambda=-1; n=100; h=1/n; c=lambda*h

-0.01

>y=sequence("(2c+5)*x[n-1]+(4c-4)*x[n-2]",[1,exp(-lambda*h)],n+1)[-1]

-2.03142847388e+057