7.1. An experiment#
Example 7.1
Solve the following initial value problem
\[
y' = -\lambda y + \lambda t + 1, \quad t\in [0, 1], ~ y(0)=1, ~\lambda=10
\]
by using the 4-step Adams-Bashforth method
\[
y_{j+1} = y_j + \frac{h}{24}\left(\,55f_j - 59f_{j-1} + 37f_{j-2} - 9f_{j-3}\right).
\]
You should
Find the accuracy of the method;
Assess the zero-stability of the method;
Examine the convergence of the method;
Implement the method in Matlab/Python, and compare the numerical solution with the exact solution
\[y=e^{-\lambda t}+t, \quad \lambda=10.\]
7.1.1. Accuracy#
We can find the error constant of the method (see Chapter 3). Rearranging the formula as
\[\begin{split}
& y_{j+1} - y_j + 0 y_{j-1} + 0 y_{j-2} + 0 y_{j-3} \\
= & h \left( 0 f_{j+1}+ \tfrac{5}{24}f_j - \tfrac{59}{24}f_{j-1} + \tfrac{37}{24}f_{j-2} - \tfrac{9}{24}f_{j-3} \right),
\end{split}\]
so
\[\begin{split}
& \alpha_0 = 0, ~
\alpha_1 = 0, ~
\alpha_2 = 0, ~
\alpha_3 = -1, ~
\alpha_4 = 1, \\
& \beta_0 = -\tfrac{9}{24}, ~
\beta_1 = \tfrac{37}{24}, ~
\beta_2= -\tfrac{59}{24}, ~
\beta_3= \tfrac{5}{24}, ~
\beta_4=0
\end{split}\]
Using Python to calculate the error constant
import math
from fractions import Fraction
a0=Fraction(0, 1);
a1=Fraction(0, 1);
a2=Fraction(0, 1);
a3=Fraction(-1, 1);
a4=Fraction(1, 1);
b0=Fraction(-9, 24);
b1=Fraction(37, 24);
b2=Fraction(-59, 24);
b3=Fraction(55, 24);
b4=Fraction(0, 1);
c0= a0+a1+a2+a3+a4;
c1= (a1+2*a2+3*a3+4*a4)-(b0+b1+b2+b3+b4);
c2= (a1+2**2*a2+3**2*a3+4**2*a4)/math.factorial(2)-(b1+2*b2+3*b3+4*b4)/math.factorial(1);
c3= (a1+2**3*a2+3**3*a3+4**3*a4)/math.factorial(3)-(b1+2**2*b2+3**2*b3+4**2*b4)/math.factorial(2);
c4= (a1+2**4*a2+3**4*a3+4**4*a4)/math.factorial(4)-(b1+2**3*b2+3**3*b3+4**3*b4)/math.factorial(3);
c5= (a1+2**5*a2+3**5*a3+4**5*a4)/math.factorial(5)-(b1+2**4*b2+3**4*b3+4**4*b4)/math.factorial(4);
print("c0=",c0)
print("c1=",c1)
print("c2=",c2)
print("c3=",c3)
print("c4=",c4)
print("c5=",c5)
c0= 0
c1= 0
c2= 0
c3= 0
c4= 0
c5= 251/720
therefore the order of accuracy of the method is \(4\), and so the method is consistent.
7.1.2. Zero-stability#
The first characteristic polynomial of the method is
\[
\rho(z) = \sum_{i=0}^{k} \alpha_i z^i = -z^3 + z^4 = -z^3 (1-z),
\]
and it has roots
\[
z_{1,2,3}=0, z_4=1.
\]
The method satisfies the root condition, so it is zero-stable.
7.1.3. Convergence#
The method is both consistent and zero-stable, so it is convergent.
7.1.4. Code Implementation#
First, let’s try a step size \(h=0.05\)
Show code cell source
import math
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['text.usetex'] = True
def f(t,lam,y):
return lam*(-y+t)+1;
def yexact(t,lam):
return math.exp(-lam*t)+t;
#interval of definition
L=1.0;
nsteps=21;
h=L/(nsteps-1);
#lambda
lam=10.0;
#declare arrays
t=np.zeros(nsteps);
y=np.zeros(nsteps);
F=np.zeros(nsteps);
yex=np.zeros(nsteps);
abserror=np.zeros(nsteps);
#initial values
t[0]=0.0;
y[0]=yexact(t[0],lam);
F[0]=f(t[0],lam,y[0]);
#calculate the first 3 starting values using the exact solution
for i in range(1,4):
t[i]=t[i-1]+h;
y[i]=yexact(t[i],lam);
F[i]=f(t[i],lam,y[i]);
#use multiple-step method to calculate y
for i in range(4,nsteps):
y[i]=y[i-1]+(h/24)*(55*F[i-1]-59*F[i-2]+37*F[i-3]-9*F[i-4]);
t[i]=t[i-1]+h;
F[i]=f(t[i],lam,y[i]);
#screen output
print(" i t y F y_ex AbsError");
for i in range(nsteps):
yex[i]=yexact(t[i],lam);
abserror[i]=abs(yex[i]-y[i]);
print("%3i %6.3f %11.5e %11.5e %11.5e %11.5e" %(i, t[i],y[i],F[i],yex[i],abserror[i]))
#exact solution for plot
ne=101;
he=L/(ne-1);
te=np.zeros(ne);
ye=np.zeros(ne);
for i in range(ne):
te[i]=i*he;
ye[i]=yexact(te[i],lam);
#plot solution
plt.figure(1)
plt.plot(t,y,'r-o',label="numerical")
plt.plot(te,ye,'b-',label="exact")
plt.xlabel("$t$", fontsize=18)
plt.ylabel("$y$", fontsize=18)
plt.legend(loc="upper center", fontsize=14)
plt.title("Solution to $y'=-10y+t+1, y(0)=1, h="+str(h)+"$",fontsize=14)
plt.yticks(fontsize=12)
plt.xticks(fontsize=12)
plt.ylim([0.2, 1.5])
plt.show()
i t y F y_ex AbsError
0 0.000 1.00000e+00 -9.00000e+00 1.00000e+00 0.00000e+00
1 0.050 6.56531e-01 -5.06531e+00 6.56531e-01 0.00000e+00
2 0.100 4.67879e-01 -2.67879e+00 4.67879e-01 0.00000e+00
3 0.150 3.73130e-01 -1.23130e+00 3.73130e-01 0.00000e+00
4 0.200 3.39611e-01 -3.96113e-01 3.35335e-01 4.27600e-03
5 0.250 3.34055e-01 1.59451e-01 3.32085e-01 1.96994e-03
6 0.300 3.56329e-01 4.36712e-01 3.49787e-01 6.54168e-03
7 0.350 3.79323e-01 7.06772e-01 3.80197e-01 8.74587e-04
8 0.400 4.26346e-01 7.36541e-01 4.18316e-01 8.03031e-03
9 0.450 4.54541e-01 9.54593e-01 4.61109e-01 6.56829e-03
10 0.500 5.19680e-01 8.03200e-01 5.06738e-01 1.29421e-02
11 0.550 5.37901e-01 1.12099e+00 5.54087e-01 1.61858e-02
12 0.600 6.27394e-01 7.26058e-01 6.02479e-01 2.49155e-02
13 0.650 6.16815e-01 1.33185e+00 6.51503e-01 3.46888e-02
14 0.700 7.51528e-01 4.84721e-01 7.00912e-01 5.06160e-02
15 0.750 6.78310e-01 1.71690e+00 7.50553e-01 7.22428e-02
16 0.800 9.04508e-01 -4.50795e-02 8.00335e-01 1.04172e-01
17 0.850 7.00699e-01 2.49301e+00 8.50203e-01 1.49505e-01
18 0.900 1.11515e+00 -1.15153e+00 9.00123e-01 2.15030e-01
19 0.950 6.41108e-01 4.08892e+00 9.50075e-01 3.08967e-01
20 1.000 1.44419e+00 -3.44187e+00 1.00005e+00 4.44142e-01
It seems the errors are too large, so let’s try a smaller step size \(h=0.025\).
Show code cell source
import math
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['text.usetex'] = True
def f(t,lam,y):
return lam*(-y+t)+1;
def yexact(t,lam):
return math.exp(-lam*t)+t;
#interval of definition
L=1.0;
nsteps=41;
h=L/(nsteps-1);
#lambda
lam=10.0;
#declare arrays
t=np.zeros(nsteps);
y=np.zeros(nsteps);
F=np.zeros(nsteps);
yex=np.zeros(nsteps);
abserror=np.zeros(nsteps);
#initial values
t[0]=0.0;
y[0]=yexact(t[0],lam);
F[0]=f(t[0],lam,y[0]);
#calculate the first 3 starting values using the exact solution
for i in range(1,4):
t[i]=t[i-1]+h;
y[i]=yexact(t[i],lam);
F[i]=f(t[i],lam,y[i]);
#use multiple-step method to calculate y
for i in range(4,nsteps):
y[i]=y[i-1]+(h/24)*(55*F[i-1]-59*F[i-2]+37*F[i-3]-9*F[i-4]);
t[i]=t[i-1]+h;
F[i]=f(t[i],lam,y[i]);
#screen output
print(" i t y F y_ex AbsError");
for i in range(nsteps):
yex[i]=yexact(t[i],lam);
abserror[i]=abs(yex[i]-y[i]);
print("%3i %6.3f %11.5e %11.5e %11.5e %11.5e" %(i, t[i],y[i],F[i],yex[i],abserror[i]))
#exact solution for plot
ne=101;
he=L/(ne-1);
te=np.zeros(ne);
ye=np.zeros(ne);
for i in range(ne):
te[i]=i*he;
ye[i]=yexact(te[i],lam);
#plot solution
plt.figure(1)
plt.plot(t,y,'r-o',label="numerical")
plt.plot(te,ye,'b-',label="exact")
plt.xlabel("$t$", fontsize=18)
plt.ylabel("$y$", fontsize=18)
plt.legend(loc="upper center", fontsize=14)
plt.title("Solution to $y'=-10y+t+1, y(0)=1, h="+str(h)+"$",fontsize=14)
plt.yticks(fontsize=12)
plt.xticks(fontsize=12)
plt.ylim([0.2, 1.5])
plt.show()
i t y F y_ex AbsError
0 0.000 1.00000e+00 -9.00000e+00 1.00000e+00 0.00000e+00
1 0.025 8.03801e-01 -6.78801e+00 8.03801e-01 0.00000e+00
2 0.050 6.56531e-01 -5.06531e+00 6.56531e-01 0.00000e+00
3 0.075 5.47367e-01 -3.72367e+00 5.47367e-01 0.00000e+00
4 0.100 4.68091e-01 -2.68091e+00 4.67879e-01 2.11274e-04
5 0.125 4.11760e-01 -1.86760e+00 4.11505e-01 2.54771e-04
6 0.150 3.73497e-01 -1.23497e+00 3.73130e-01 3.66798e-04
7 0.175 3.49106e-01 -7.41055e-01 3.48774e-01 3.31602e-04
8 0.200 3.35702e-01 -3.57017e-01 3.35335e-01 3.66386e-04
9 0.225 3.30703e-01 -5.70254e-02 3.30399e-01 3.03320e-04
10 0.250 3.32393e-01 1.76066e-01 3.32085e-01 3.08442e-04
11 0.275 3.39173e-01 3.58274e-01 3.38928e-01 2.44736e-04
12 0.300 3.50027e-01 4.99728e-01 3.49787e-01 2.40123e-04
13 0.325 3.63959e-01 6.10410e-01 3.63774e-01 1.84789e-04
14 0.350 3.80376e-01 6.96242e-01 3.80197e-01 1.78429e-04
15 0.375 3.98651e-01 7.63486e-01 3.98518e-01 1.33675e-04
16 0.400 4.18444e-01 8.15558e-01 4.18316e-01 1.28560e-04
17 0.425 4.39358e-01 8.56420e-01 4.39264e-01 9.38065e-05
18 0.450 4.61200e-01 8.88003e-01 4.61109e-01 9.06607e-05
19 0.475 4.83716e-01 9.12840e-01 4.83652e-01 6.43233e-05
20 0.500 5.06801e-01 9.31991e-01 5.06738e-01 6.29575e-05
21 0.525 5.30291e-01 9.47092e-01 5.30248e-01 4.32860e-05
22 0.550 5.54130e-01 9.58700e-01 5.54087e-01 4.32345e-05
23 0.575 5.78211e-01 9.67886e-01 5.78183e-01 2.86609e-05
24 0.600 6.02508e-01 9.74918e-01 6.02479e-01 2.94545e-05
25 0.625 6.26949e-01 9.80508e-01 6.26930e-01 1.86975e-05
26 0.650 6.51523e-01 9.84766e-01 6.51503e-01 1.99579e-05
27 0.675 6.76183e-01 9.88171e-01 6.76171e-01 1.20220e-05
28 0.700 7.00925e-01 9.90746e-01 7.00912e-01 1.34789e-05
29 0.725 7.25718e-01 9.92822e-01 7.25710e-01 7.61376e-06
30 0.750 7.50562e-01 9.94378e-01 7.50553e-01 9.09084e-06
31 0.775 7.75435e-01 9.95645e-01 7.75431e-01 4.74128e-06
32 0.800 8.00342e-01 9.96584e-01 8.00335e-01 6.13383e-06
33 0.825 8.25264e-01 9.97358e-01 8.25261e-01 2.89364e-06
34 0.850 8.50208e-01 9.97924e-01 8.50203e-01 4.14733e-06
35 0.875 8.75160e-01 9.98398e-01 8.75158e-01 1.72106e-06
36 0.900 9.00126e-01 9.98738e-01 9.00123e-01 2.81458e-06
37 0.925 9.25097e-01 9.99029e-01 9.25096e-01 9.87817e-07
38 0.950 9.50077e-01 9.99232e-01 9.50075e-01 1.92015e-06
39 0.975 9.75059e-01 9.99412e-01 9.75058e-01 5.37202e-07
40 1.000 1.00005e+00 9.99533e-01 1.00005e+00 1.31874e-06
The solution looks quite good now.
Questions
When the step size is \(h=0.05\), the numerical solution is unstable. After we reduced it to \(h=0.0025\), the solution becomes stable.
Shall we limit the step size to keep the solution stable?
What is the maximum step size we can use?
If we use a different linear multistep method, will the maximum allowable step size be the same?