22======= Analysis of schemes for the diffusion equation =======
33label{diffu:pde1:analysis}
44
5+ idx{noise! sawtooth-like}
6+
57The numerical experiments in Sections
68ref{diffu:pde1:FE:experiments} and ref{diffu:pde1:theta:experiments}
79reveal that there are some
@@ -25,6 +27,9 @@ Figures ref{diffu:pde1:FE:fig:F=0.5}-ref{diffu:pde1:CN:fig:F=10}.
2527===== Properties of the solution =====
2628label{diffu:pde1:analysis:uex}
2729
30+ idx{error function (erf)!} idx{step function}
31+ idx{error function (erf)! complementary}
32+
2833A particular characteristic of diffusive processes, governed
2934by an equation like
3035
@@ -102,11 +107,13 @@ label{diffu:analysis:pde1:p1:erf:uR}
102107\tp
103108\end{align}
104109!et
105- For small enough $t$, $u(0,t)\approx 1$ and $u(1,t)\approx 1 $, but as
110+ For small enough $t$, $u(0,t)\approx 1$ and $u(1,t)\approx 0 $, but as
106111$t\rightarrow\infty$, $u(x,t)\rightarrow 1/2$ on $[0,1]$.
107112
108113=== Solution for a Gaussian pulse ===
109114
115+ idx{Dirac delta function}
116+
110117The standard diffusion equation $u_t = \dfc u_{xx}$ admits a
111118Gaussian function as solution:
112119
@@ -132,6 +139,8 @@ the pulse diffuse and flatten out.
132139
133140=== Solution for a sine component ===
134141
142+ idx{smoothing} idx{noise! removing} idx{signal processing}
143+
135144Also, (ref{diffu:pde1:eq}) admits a solution of the form
136145
137146!bt
@@ -150,8 +159,7 @@ inserting (ref{diffu:pde1:sol1}) in (ref{diffu:pde1:eq}) gives the constraint
150159\end{equation*}
151160!et
152161
153-
154- A very important feature is that the initial shape $I(x)=Q\sin kx$
162+ A very important feature is that the initial shape $I(x)=Q\sin\left( kx\right)$
155163undergoes a damping $\exp{(-\dfc k^2t)}$, meaning that
156164rapid oscillations in space, corresponding to large $k$, are very much
157165faster dampened than slow oscillations in space, corresponding to small
@@ -273,6 +281,8 @@ denotes $u$ at time $t_n$.
273281
274282=== Stability ===
275283
284+ idx{stability}
285+
276286The exact amplification factor is $\Aex=\exp{(-\dfc^2 k^2\Delta t)}$.
277287We should therefore require $|A| < 1$ to have a decaying numerical
278288solution as well. If
@@ -284,6 +294,8 @@ idx{amplification factor}
284294
285295=== Accuracy ===
286296
297+ idx{accuracy}
298+
287299To determine how accurately a finite difference scheme treats one
288300wave component (ref{diffu:pde1:analysis:uni}), we see that the basic
289301deviation from the exact solution is reflected in how well
@@ -297,6 +309,8 @@ make Taylor expansions of $A/\Aex$ to see the error more analytically.
297309
298310=== Truncation error ===
299311
312+ idx{diffusion equation! truncation error}
313+
300314As an alternative to examining the accuracy of the damping of a wave
301315component, we can perform a general truncation error analysis as
302316explained in ref[Appendix ref{ch:trunc}][ in
@@ -312,6 +326,8 @@ verifying codes based on empirical estimation of convergence rates.
312326===== Analysis of the Forward Euler scheme =====
313327label{diffu:pde1:analysis:FE}
314328
329+ idx{diffusion equation! numerical Fourier number}
330+
315331# 2DO: refer to vib and wave
316332
317333
@@ -402,6 +418,8 @@ The method hence becomes very expensive for fine spatial meshes.
402418
403419=== Accuracy ===
404420
421+ idx{Taylor series}
422+
405423Since $A$ is expressed in terms of $F$ and the parameter we now call
406424$p=k\Delta x/2$, we should also express $\Aex$ by $F$ and $p$. The exponent
407425in $\Aex$ is $-\dfc k^2\Delta t$, which equals $-F k^2\Delta x^2=-F4p^2$.
@@ -578,7 +596,7 @@ The exact numerical solution is hence
578596
579597!bt
580598\begin{equation}
581- u^n_q = \left(\frac{ 1 - 2F\sin^2p}{1 + 2F\sin^2p}\right)^ne^{ikp \Delta x}
599+ u^n_q = \left(\frac{ 1 - 2F\sin^2p}{1 + 2F\sin^2p}\right)^ne^{ikq \Delta x}
582600\tp
583601\end{equation}
584602!et
@@ -605,21 +623,23 @@ cite{Langtangen_deqbook_trunc}]:
605623===== Analysis of the Leapfrog scheme =====
606624label{diffu:pde1:analysis:leapfrog}
607625
626+ idx{Leapfrog scheme}
627+
608628An attractive feature of the Forward Euler scheme is the explicit
609629time stepping and no need for solving linear systems. However, the
610630accuracy in time is only $\Oof{\Delta t}$. We can get an explicit
611631*second-order* scheme in time by using the Leapfrog method:
612632
613633!bt
614- \[ [D_{2t} u = \dfc D_xDx u + f]^n_i \tp\]
634+ \[ [D_{2t} u = \dfc D_xDx u + f]^n_q \tp\]
615635!et
616636Written out,
617637
618638!bt
619- \[ u ^{n+1} = u ^{n-1} + \frac{2\dfc\Delta t}{\Delta x^2}
620- (u^{n}_{i +1} - 2u^n_i + u^n_{i -1}) + f(x_i ,t_n)\tp\]
639+ \[ u_q ^{n+1} = u_q ^{n-1} + \frac{2\dfc\Delta t}{\Delta x^2}
640+ (u^{n}_{q +1} - 2u^n_q + u^n_{q -1}) + f(x_q ,t_n)\tp\]
621641!et
622- We need some formula for the first step, $u^1_i $, but for that we can use
642+ We need some formula for the first step, $u^1_q $, but for that we can use
623643a Forward Euler step.
624644
625645Unfortunately, the Leapfrog scheme is always unstable for the
@@ -690,7 +710,8 @@ respectively, and we see how short waves pollute the overall solution.
690710===== Analysis of the 2D diffusion equation =====
691711label{diffu:2D:analysis}
692712
693- We first consider the 2D diffusion equation
713+ Diffusion in several dimensions is treated later, but it is appropriate to
714+ include the analysis here. We first consider the 2D diffusion equation
694715
695716!bt
696717\[ u_{t} = \dfc(u_{xx} + u_{yy}),\]
@@ -711,7 +732,7 @@ and the schemes have discrete versions of this Fourier component:
711732For the Forward Euler discretization,
712733
713734!bt
714- \[ [D_t^+u = \dfc(D_xD_x u + D_yD_y u)]_{i,j }^n,\]
735+ \[ [D_t^+u = \dfc(D_xD_x u + D_yD_y u)]_{q,r }^n,\]
715736!et
716737we get
717738
@@ -752,7 +773,7 @@ The complete numerical solution for a wave component is
752773!bt
753774\begin{equation}
754775u^{n}_{q,r} = A(1 - 4F_x\sin^2 p_x - 4F_y\sin^2 p_y)^n
755- e^{i(k_xp \Delta x + k_yq \Delta y)}\tp
776+ e^{i(k_xq \Delta x + k_yr \Delta y)}\tp
756777label{diffu:2D:analysis:FE:numexact}
757778\end{equation}
758779!et
@@ -781,7 +802,7 @@ where $d$ is the number of space dimensions: $d=1,2,3$.
781802The Backward Euler method,
782803
783804!bt
784- \[ [D_t^-u = \dfc(D_xD_x u + D_yD_y u)]_{i,j }^n,\]
805+ \[ [D_t^-u = \dfc(D_xD_x u + D_yD_y u)]_{q,r }^n,\]
785806!et
786807results in
787808
@@ -809,9 +830,9 @@ label{diffu:2D:analysis:BN:numexact}
809830With a Crank-Nicolson discretization,
810831
811832!bt
812- \[ [D_tu]^{n+\half}_{i,j } =
813- \half [\dfc(D_xD_x u + D_yD_y u)]_{i,j }^{n+1} +
814- \half [\dfc(D_xD_x u + D_yD_y u)]_{i,j }^n,\]
833+ \[ [D_tu]^{n+\half}_{q,r } =
834+ \half [\dfc(D_xD_x u + D_yD_y u)]_{q,r }^{n+1} +
835+ \half [\dfc(D_xD_x u + D_yD_y u)]_{q,r }^n,\]
815836!et
816837we have, after some algebra,
817838
@@ -853,8 +874,8 @@ label{diffu:2D:analysis:CN:numexact}
853874
854875The behavior of the solution generated by Forward Euler discretization in time (and centered
855876differences in space) is summarized at the end of
856- Section ref{diffu:pde1:FE:experiments}. Can we from the analysis
857- above explain the behavior?
877+ Section ref{diffu:pde1:FE:experiments}. Can we, from the analysis
878+ above, explain the behavior?
858879
859880We may start by looking at Figure ref{diffu:pde1:FE:fig:F=0.51}
860881where $F=0.51$. The figure shows that the solution is unstable and
0 commit comments