Instability in a multistep method¶

Consider the ridiculously simple IVP \(u'=u\), \(u(0)=1\), whose solution is \(e^t\).

using FundamentalsNumericalComputation

dudt = (u,t) -> u;
u_exact = exp;
a = 0.0;  b = 1.0;

Let’s apply the LIAF method to this problem for varying fixed step sizes. We’ll measure the error at the time \(t=1\).

n = [5,10,20,40,60]
err = zeros(size(n))
t = [];  u = [];
for (j,n) = enumerate(n)
    h = (b-a)/n
    t = [ a + i*h for i=0:n ]
    u = [1; u_exact(h); zeros(n-1)];
    f = [dudt(u[1],t[1]); zeros(n)];
    for i = 2:n
        f[i] = dudt(u[i],t[i])
        u[i+1] = -4*u[i] + 5*u[i-1] + h*(4*f[i]+2*f[i-1])
    end
    err[j] = abs(u_exact(b) - u[end])
end

pretty_table( (n=n,h=(b-a)./n,error=err), backend=:html )

n	h	error
Int64	Float64	Float64
5	0.2	0.01604519922158687
10	0.1	2.845479238377764
20	0.05	1.6224961462218664e6
40	0.025	9.3442027524627e18
60	0.016666666666666666	1.740134731635269e32

The error starts out promisingly, but things explode from there. A graph of the last numerical attempt yields a clue.

plot(t,abs.(u),m=:o,label="",
    xlabel="t",yaxis=(:log10,"|u|"),title="LIAF solution")

It’s clear that the solution is growing exponentially in time.

Fundamentals of Numerical Computation

Instability in a multistep method¶