MATH761

MATH76109/07/2018

Lab03: Analysis of finite difference schemes

1Finite difference approximations

Mathematica implementation of finite difference calculus

In[1]:=

x \[Element] Reals; $Assumptions = h > 0;

Dp[f_[x_]] := f[x + h] - f[x];

Dm[f_[x_]] := f[x] - f[x - h];

d[f_[x_]] := f[x+h/2]-f[x-h/2];

Dp[a_ f_[x_]]=a Dp[f[x]]; Dp[a_+b_]=Dp[a]+Dp[b]; Dp[a_ (b_+c_)]=a Dp[b]+a Dp[c]; Dp[a_]=0;

Dm[a_ f_[x_]]=a Dm[f[x]]; Dm[a_+b_]=Dm[a]+Dm[b]; Dm[a_ (b_+c_)]=a Dm[b]+a Dm[c]; Dm[a_]=0;

d[a_ f_[x_]]=a d[f[x]]; d[a_+b_]= d[a]+ d[b]; d[a_ (b_+c_)]=a d[b]+a d[c]; d[a_]=0;

fdSort[expr_]:=Assuming[{h > 0, x > 0}, Sort[Level[expr, 1], Level[#1, {Depth[#1] - 1}] < Level[#2, {Depth[#2] - 1}] &]];

Q[x_]=Subscript[Q,x];

Dp[f_[x_], p_] :=

Together[Normal[1/h Series[Log[1 + z], {z, 0, p}]] /.

Table[z^j -> Nest[Dp, f[x], j], {j, 1, p}]];

Dm[f_[x_], p_] :=

Together[Normal[-1/h Series[Log[1 - z], {z, 0, p}]] /.

Table[z^j -> Nest[Dm, f[x], j], {j, 1, p}]];

d[f_[x_], p_] := Together[Normal[2/h Series[ArcSinh[z/2], {z, 0, 2 p - 1}]] /.

Table[z^j -> Nest[d, f[x], j], {j, 1, 2 p - 1}]];

Dp[f_[x_], p_, k_] :=

Together[Expand[(Normal[1/h Series[Log[1 + z], {z, 0, p}]])^k] /.

Table[z^j -> Nest[Dp, f[x], j], {j, 1, p k}]];

Dm[f_[x_], p_, k_] :=

Together[Expand[(Normal[-1/h Series[Log[1 - z], {z, 0, p}]])^k] /.

Table[z^j -> Nest[Dm, f[x], j], {j, 1, p k}]];

d[f_[x_], p_, k_] := Together[Expand[(Normal[2/h Series[ArcSinh[z/2], {z, 0, 2 p - 1}]])^k]

/. Table[z^j -> Nest[d, f[x], j], {j, 1, (2 p - 1)k}]];

{qt1,qxx2,qxx4}={Dp[f[x],1],d[f[x],1,2],d[f[x],2,2]}

In[2]:=

2Interpolating polynomial

The Newton form of the interpolating polynomial of data $𝒟 = {(t^{n} = n h, Q^{n}), n = 0, 1, \dots, r}$ is

p_{r} (t) = Q^{0} + [Q^{1}, Q^{0}] (t - t^{0}) + \dots + [Q^{r}, \dots, Q^{1}, Q^{0}] (t - t^{0}) \dots (t - t^{r - 1})

In[2]:=

pNewton[q_[t_],r_] := Sum[ Nest[Dp,q[0],j]/h^j Product[(t - n h),{n,0,j-1}]/j!,{j,0,r}];

Table[Simplify[pNewton[q[t],r]],{r,0,3}] // TableForm

$\begin{array}{c} q (0) \\ \frac{t (q (h) - q (0))}{h} + q (0) \\ \frac{t (- 2 q (h) + q (2 h) + q (0)) (t - h)}{2 h^{2}} + \frac{t (q (h) - q (0))}{h} + q (0) \\ \frac{3 t q (h) (6 h^{2} - 5 h t + t^{2}) + (h - t) (3 t q (2 h) (t - 3 h) + (2 h - t) (t q (3 h) + 3 h q (0) - q (0) t))}{6 h^{3}} \end{array}$

In[3]:=

Verify the above implementation for the first few interpolating polynomials

In[3]:=

p1[t_] = pNewton[q[t],1]

$\frac{t (q (h) - q (0))}{h} + q (0)$

In[4]:=

{p1[0],p1[h]}

${q (0), q (h)}$

In[5]:=

p2[t_] = pNewton[q[t],2]

$\frac{t (- 2 q (h) + q (2 h) + q (0)) (t - h)}{2 h^{2}} + \frac{t (q (h) - q (0))}{h} + q (0)$

In[6]:=

Simplify[{p2[0],p2[h],p2[2h]}]

${q (0), q (h), q (2 h)}$

In[7]:=

p3[t_] = pNewton[q[t],3]

$\frac{t (q (h) - q (2 h) - 2 (q (2 h) - q (h)) + q (3 h) - q (0)) (t - 2 h) (t - h)}{6 h^{3}} + \frac{t (- 2 q (h) + q (2 h) + q (0)) (t - h)}{2 h^{2}} + \frac{t (q (h) - q (0))}{h} + q (0)$

In[8]:=

Simplify[{p3[0],p3[h],p3[2h],p3[3h]}]

${q (0), q (h), q (2 h), q (3 h)}$

In[9]:=

3Linear multistep methods

In a linear multistep method (LMM) the solution to ODE $q^{'} = f (q, t)$ is sought by linear combination of values $Q^{n} ≅ q (t^{n})$ , $F^{n} ≅ (t^{n}, q (t^{n}))$ with $t^{n} = n k$ .

3.1Adams-Bashforth

Adams-Bashforth schemes are obtained by integrating the ODE over interval $[t^{n + r}, t^{n + r - 1}]$

Q^{n + 1 + r} = Q^{n + r} + \int_{(n + r) k}^{(n + 1 + r) k} f (t, q (t)) d t,

and approximating $f$ by the interpolating polynomial $p_{r} (t)$ . The following computes the lhs, rhs of the above formula, setting $n = 0$ .

In[9]:=

AdamsBashforth[q_[t_],r_] :=

{q[(r+1)k],q[r k]+Integrate[ pNewton[f[t],r] /. h->k, {t,r k,(r+1) k}]};

ABschemes = Table[AdamsBashforth[q[t],r],{r,0,4}] // Apart;

ABschemes // TableForm

$\begin{array}{cc} q (k) & f (0) k + q (0) \\ q (2 k) & q (k) - \frac{1}{2} k (f (0) - 3 f (k)) \\ q (3 k) & \frac{1}{12} k (- 16 f (k) + 23 f (2 k) + 5 f (0)) + q (2 k) \\ q (4 k) & q (3 k) - \frac{1}{24} k (- 37 f (k) + 59 f (2 k) - 55 f (3 k) + 9 f (0)) \\ q (5 k) & \frac{1}{720} k (- 1274 f (k) + 2616 f (2 k) - 2774 f (3 k) + 1901 f (4 k) + 251 f (0)) + q (4 k) \end{array}$

In[10]:=

3.2Adams-Moulton

Adams-Moulton are are obtained by integrating the ODE over interval $[t^{n + r - 1}, t^{n + r}]$

Q^{n + r} = Q^{n + r - 1} + \int_{(n + r - 1) k}^{(n + r) k} f (t, q (t)) d t,

and approximating $f$ by the interpolating polynomial $p_{r} (t)$ . The following computes the lhs, rhs of the above formula, setting $n = 0$ .

In[10]:=

AdamsMoulton[q_[t_],r_] :=

{q[r k],q[(r-1) k]+Integrate[ pNewton[f[t],r] /. h->k, {t,(r-1) k,r k}]};

AMschemes = Table[AdamsMoulton[q[t],r],{r,0,4}] // Apart;

AMschemes // TableForm

$\begin{array}{cc} q (0) & f (0) k + q (- k) \\ q (k) & \frac{1}{2} k (f (k) + f (0)) + q (0) \\ q (2 k) & q (k) - \frac{1}{12} k (- 8 f (k) - 5 f (2 k) + f (0)) \\ q (3 k) & \frac{1}{24} k (- 5 f (k) + 19 f (2 k) + 9 f (3 k) + f (0)) + q (2 k) \\ q (4 k) & q (3 k) - \frac{1}{720} k (- 106 f (k) + 264 f (2 k) - 646 f (3 k) - 251 f (4 k) + 19 f (0)) \end{array}$

In[11]:=

4Stability analysis

Define a function to obtain the characteristic polynomial given as a pair {lhs,rhs}

In[11]:=

CharPoly[s_] := Module[{cp,g,F,Q},

g[t_]=f[q[t]]; F[q_] = lambda q; Q[t_] = zeta^(t/k);

cp = Simplify[s[[1]] - s[[2]] /. f->g /. f->F /. q->Q];

Return[Simplify[cp /. k->(z/lambda)]]

];

CharPoly[ABschemes[[1]]]

$- z + ζ - 1$

In[12]:=

Compute and display the characteristic polynomials for Adams-Bashforth, Adams-Moulton schemes

In[12]:=

ABCPs = Map[CharPoly[#]&,ABschemes];

ABCPs // TableForm

$\begin{array}{c} - z + ζ - 1 \\ (ζ - 1) ζ - \frac{1}{2} z (3 ζ - 1) \\ (ζ - 1) ζ^{2} - \frac{1}{12} z (23 ζ^{2} - 16 ζ + 5) \\ (ζ - 1) ζ^{3} - \frac{1}{24} z (55 ζ^{3} - 59 ζ^{2} + 37 ζ - 9) \\ (ζ - 1) ζ^{4} - \frac{1}{720} z (1901 ζ^{4} - 2774 ζ^{3} + 2616 ζ^{2} - 1274 ζ + 251) \end{array}$

In[13]:=

AMCPs = Map[CharPoly[#]&,AMschemes];

AMCPs // TableForm

$\begin{array}{c} - z - \frac{1}{ζ} + 1 \\ - \frac{1}{2} z (ζ + 1) + ζ - 1 \\ (ζ - 1) ζ - \frac{1}{12} z (5 ζ^{2} + 8 ζ - 1) \\ (ζ - 1) ζ^{2} - \frac{1}{24} z (9 ζ^{3} + 19 ζ^{2} - 5 ζ + 1) \\ (ζ - 1) ζ^{3} - \frac{1}{720} z (251 ζ^{4} + 646 ζ^{3} - 264 ζ^{2} + 106 ζ - 19) \end{array}$

In[14]:=

Define a function to extract the boundary locus $| ζ | = 1$ , i.e. $ζ = e^{i θ}$

In[14]:=

BLocus[p_] := z /. Solve[ p == 0, z][[1]] /. zeta->Exp[I theta];

BLocus[ABCPs[[1]]]

$- 1 + e^{i θ}$

In[15]:=

Compute and display the boundary locii for Adams-Bashforth, Adams-Moulton schemes

In[15]:=

ABBLs = Map[BLocus[#]&, ABCPs];

ABBLs // TableForm

$\begin{array}{c} - 1 + e^{i θ} \\ \frac{2 e^{i θ} (- 1 + e^{i θ})}{- 1 + 3 e^{i θ}} \\ \frac{12 e^{2 i θ} (- 1 + e^{i θ})}{- 16 e^{i θ} + 23 e^{2 i θ} + 5} \\ \frac{24 e^{3 i θ} (- 1 + e^{i θ})}{37 e^{i θ} - 59 e^{2 i θ} + 55 e^{3 i θ} - 9} \\ \frac{720 e^{4 i θ} (- 1 + e^{i θ})}{- 1274 e^{i θ} + 2616 e^{2 i θ} - 2774 e^{3 i θ} + 1901 e^{4 i θ} + 251} \end{array}$

In[16]:=

AMBLs = Map[BLocus[#]&, AMCPs];

AMBLs // TableForm

$\begin{array}{c} e^{- i θ} (- 1 + e^{i θ}) \\ \frac{2 (- 1 + e^{i θ})}{1 + e^{i θ}} \\ \frac{12 e^{i θ} (- 1 + e^{i θ})}{8 e^{i θ} + 5 e^{2 i θ} - 1} \\ \frac{24 e^{2 i θ} (- 1 + e^{i θ})}{- 5 e^{i θ} + 19 e^{2 i θ} + 9 e^{3 i θ} + 1} \\ \frac{720 e^{3 i θ} (- 1 + e^{i θ})}{106 e^{i θ} - 264 e^{2 i θ} + 646 e^{3 i θ} + 251 e^{4 i θ} - 19} \end{array}$

Define a function to plot a family of boundary locii, and apply it to Adams-Bashforth and Adams-Moulton schemes

In[17]:=

BLocusPlot[bl_,blname_] :=

ParametricPlot[Evaluate[Map[{Re[#],Im[#]}&,bl]],{theta,0,2Pi},Axes->True,Frame->True,

FrameLabel->{"Re(z)","Im(z)"},GridLines->Automatic,AspectRatio->Automatic,

PlotLabel->blname,PlotLegends->Automatic,PlotRange->{{-3,3},{-3,3}}];

ABBLplots = BLocusPlot[ABBLs,"Adams-Bashforth boundary locii"];

Export["/home/student/courses/MATH761/AdamsBashforthBoundaryLocii.pdf",ABBLplots]

$/home/student/courses/MATH761/AdamsBashforthBoundaryLocii.pdf$

In[18]:=

AMBLplots = BLocusPlot[AMBLs,"Adams-Moulton boundary locii"];

Export["/home/student/courses/MATH761/AdamsMoultonBoundaryLocii.pdf",AMBLplots]

$/home/student/courses/MATH761/AdamsMoultonBoundaryLocii.pdf$

In[19]:=

Figure 1. Boundary locii for LMMs