Stability of BDF


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

>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
>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
>function ystep(y0,y1,yd2,h) &= ps(h)|ratsimp
2 h yd2 + 4 y1 - y0
-------------------
3

>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
>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

>difgl &= solve(ystep(y0,y1,lambda*y2,h)=y2,y2)
y0 - 4 y1
[y2 = --------------]
2 h lambda - 3
>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
>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


>function q(z) &= 3-4*z+z^2
2
z - 4 z + 3
>z=exp(I*linspace(0,2pi,1000)); plot2d(q(z),r=10); insimg;

>&realpart(q(cos(x)+I*sin(x)))|trigsimp
2 2
- sin (x) + cos (x) - 4 cos(x) + 3
>plot2d("2*x^2-4*x+2",-1,1); insimg;

>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

>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
>plot2d(q(z),r=50); insimg;




>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
>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
>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
>y2c &= ratsimp(y2 with lambda=c/h)
(4 c - 4) y1 + (2 c + 5) y0
>&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
>z=exp(I*linspace(0,2pi,1000)); ...
>plot2d(h(z),r=6); insimg;

>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
Examples Homepage