Merge branch 'master' of github.com:hplgit/fdm-book

Author: sveinlin

doc/.src/chapters/nonlin/nonlin_split.do.txt

@@ -297,15 +297,16 @@ splitting methods.
 label{nonlin:splitting:RD_linearR}
 
 The methods above may be explored in detail through a specific
-computational example in which we compute the convergence
-rates associated with four different solution approaches for the
-reaction-diffusion equation with a linear reaction term, i.e. $f(u)=-bu$. 
-The methods comprise solving without splitting (just straight Forward Euler), 
-ordinary splitting, first order Strang splitting, and second order Strang splitting.
-In all four methods, a standard centered difference approximation is
-used for the spatial second derivative. The methods share the error model
-$E = C h^r$, while differing in the step $h$ (being either $dx^2$ or $dx$) and 
-the convergence rate $r$ (being either 1 or 2). 
+computational example in which we compute the convergence rates
+associated with four different solution approaches for the
+reaction-diffusion equation with a linear reaction term,
+i.e. $f(u)=-bu$.  The methods comprise solving without splitting (just
+straight Forward Euler), ordinary splitting, first order Strang
+splitting, and second order Strang splitting.  In all four methods, a
+standard centered difference approximation is used for the spatial
+second derivative. The methods share the error model $E = C h^r$,
+while differing in the step $h$ (being either $\Delta x^2$ or $\Delta
+x$) and the convergence rate $r$ (being either 1 or 2).
 
 All code commented below is found in the file
 "`split_diffu_react.py`": "${src_nonlin}/split_diffu_react.py". When executed,
@@ -314,65 +315,78 @@ rate computations are handled:
 
 @@@CODE src-nonlin/split_diffu_react.py fromto: def convergence_rates@
 
-Now, with respect to the error ($E = C h^r$), the Forward Euler scheme, the 
-ordinary splitting scheme and first order Strang splitting scheme are all
-first order ($r = 1$), with a step $h = \Delta x^2 = K^{-1}\Delta t$, where 
-$K$ is some constant. This implies that the *ratio* $\frac{\Delta t}{\Delta x^2}$ 
-must be held constant during convergence rate calculations. Furthermore, 
-the Fourier number $F = \frac{\alpha \Delta t}{\Delta x^2}$ is upwards limited to 
-$F = 0.5$, being the stability limit with explicit schemes. Thus, in these cases,
-we use the fixed value of $F$ and a given (but changing) spatial resolution $\Delta x$ to
-compute the corresponding value of $\Delta t$ according to the expression for $F$.
-This assures that $\frac{\Delta t}{\Delta x^2}$ is kept constant. The loop in 
-`convergence_rates` runs over a chosen set of grid points (`Nx_values`) which gives 
-a doubling of spatial resolution with each iteration ($\Delta x$ is halved).
-
-For the second order Strang splitting scheme, we have $r = 2$ and a step 
-$h = \Delta x = K^{-1}\Delta t$, where $K$ again is some constant. In this case,
-it is thus the ratio $\frac{\Delta t}{\Delta x}$ that must be held constant 
-during the convergence rate calculations. From the expression for $F$, it is 
-clear then that $F$ must change with each halving of $\Delta x$. In fact, if 
-$F$ is doubled each time $\Delta x$ is halved, the ratio $\frac{\Delta t}{\Delta x}$ 
-will be constant (this follows, e.g., from the expression for $F$). This us utilized
-in our code.
-
-A solver `diffusion_theta` is used in each of the four solution approaches: 
+Now, with respect to the error ($E = C h^r$), the Forward Euler
+scheme, the ordinary splitting scheme and first order Strang splitting
+scheme are all first order ($r = 1$), with a step $h = \Delta x^2 =
+K^{-1}\Delta t$, where $K$ is some constant. This implies that the
+*ratio* $\frac{\Delta t}{\Delta x^2}$ must be held constant during
+convergence rate calculations. Furthermore, the Fourier number $F =
+\frac{\alpha \Delta t}{\Delta x^2}$ is upwards limited to $F = 0.5$,
+being the stability limit with explicit schemes. Thus, in these cases,
+we use the fixed value of $F$ and a given (but changing) spatial
+resolution $\Delta x$ to compute the corresponding value of $\Delta t$
+according to the expression for $F$.  This assures that $\frac{\Delta
+t}{\Delta x^2}$ is kept constant. The loop in `convergence_rates` runs
+over a chosen set of grid points (`Nx_values`) which gives a doubling
+of spatial resolution with each iteration ($\Delta x$ is halved).
+
+For the second order Strang splitting scheme, we have $r = 2$ and a
+step $h = \Delta x = K^{-1}\Delta t$, where $K$ again is some
+constant. In this case, it is thus the ratio $\frac{\Delta t}{\Delta
+x}$ that must be held constant during the convergence rate
+calculations. From the expression for $F$, it is clear then that $F$
+must change with each halving of $\Delta x$. In fact, if $F$ is
+doubled each time $\Delta x$ is halved, the ratio $\frac{\Delta
+t}{\Delta x}$ will be constant (this follows, e.g., from the
+expression for $F$). This is utilized in our code.
+
+A solver `diffusion_theta` is used in each of the four solution approaches:
 
 @@@CODE src-nonlin/split_diffu_react.py fromto: def diffusion_theta@def reaction_FE
 
-For the no splitting approach with Forward Euler in time, this solver handles both
-the diffusion and the reaction term. When splitting, `diffusion_theta` takes care 
-of the diffusion term only, while the reaction
-term is handled either by a Forward Euler scheme in `reaction_FE`, or by a
-second order Adams-Bashforth scheme from Odespy. The `reaction_FE` function
-covers one complete time step `dt` during ordinary splitting, while Strang
-splitting (both first and second order) applies it with `dt/2` twice during 
-each time step `dt`. Since the reaction term typically represents a much faster
-process than the diffusion term, a further refinement of the time step is made
-possible in `reaction_FE`. It was implemented as
+For the no splitting approach with Forward Euler in time, this solver
+handles both the diffusion and the reaction term. When splitting,
+`diffusion_theta` takes care of the diffusion term only, while the
+reaction term is handled either by a Forward Euler scheme in
+`reaction_FE`, or by a second order Adams-Bashforth scheme from
+Odespy. The `reaction_FE` function covers one complete time step `dt`
+during ordinary splitting, while Strang splitting (both first and
+second order) applies it with `dt/2` twice during each time step
+`dt`. Since the reaction term typically represents a much faster
+process than the diffusion term, a further refinement of the time step
+is made possible in `reaction_FE`. It was implemented as
 
 @@@CODE src-nonlin/split_diffu_react.py fromto: def reaction_FE@def ordinary_splitting
 
-With the ordinary splitting approach, each time step `dt` is covered twice. First 
-computing the impact of the reaction term, then the contribution from the diffusion term:
+With the ordinary splitting approach, each time step `dt` is covered
+twice. First computing the impact of the reaction term, then the
+contribution from the diffusion term:
 
 @@@CODE src-nonlin/split_diffu_react.py fromto: def ordinary_splitting@def Strang_splitting_1st
 
-For the two Strang splitting approaches, each time step `dt` is handled by first computing
-the reaction step for (the first) `dt/2`, followed by a diffusion step `dt`, before the reaction step
-is treated once again for (the remaining) `dt/2`. Since first order Strang splitting is no better than first order accurate, both the reaction and diffusion steps are computed explicitly. The 
-solver was implemented as 
+For the two Strang splitting approaches, each time step `dt` is
+handled by first computing the reaction step for (the first) `dt/2`,
+followed by a diffusion step `dt`, before the reaction step is treated
+once again for (the remaining) `dt/2`. Since first order Strang
+splitting is no better than first order accurate, both the reaction
+and diffusion steps are computed explicitly. The solver was
+implemented as
 
 @@@CODE src-nonlin/split_diffu_react.py fromto: def Strang_splitting_1st@def Strang_splitting_2nd
 
-The second order version of the Strang splitting approach utilizes a second order Adams-Bashforth
-solver for the reaction part and a Crank-Nicolson scheme for the diffusion part. The solver has the
-same structure as the one for first order Strang splitting and was implemented as
+The second order version of the Strang splitting approach utilizes a
+second order Adams-Bashforth solver for the reaction part and a
+Crank-Nicolson scheme for the diffusion part. The solver has the same
+structure as the one for first order Strang splitting and was
+implemented as
 
 @@@CODE src-nonlin/split_diffu_react.py fromto: def Strang_splitting_2nd@def convergence_rates
 
-When executing `split_diffu_react.py`, we find that the estimated convergence rates
-are as expected. The second order Strang splitting gives the least error (about $4e^{-5}$) and has second order convergence ($r = 2$), while the remaining three approaches have first order convergence ($r = 1$).
+When executing `split_diffu_react.py`, we find that the estimated
+convergence rates are as expected. The second order Strang splitting
+gives the least error (about $4e^{-5}$) and has second order
+convergence ($r = 2$), while the remaining three approaches have first
+order convergence ($r = 1$).
 
 
 
@@ -411,7 +425,7 @@ diffusion problem solved by a Forward Euler method one has
 where $F=\dfc\Delta t/\Delta x^2$ is the mesh Fourier number and $p=k\Delta x/2$
 is a dimensionless number reflecting the spatial resolution (number of points
 per wave length in space). For the decay problem $u'=-\beta u$, we have
-$A=1 - q$, where $q$ is a dimensionless parameter reflecting the resolution
+$A=1 - q$, where $q$ is a dimensionless parameter for the resolution
 in the decay problem: $q = \beta\Delta t$.
 
 The original model problem can also be discretized by a Forward Euler scheme,
@@ -443,5 +457,5 @@ Then we use this as input to the decay algorithm and arrive at
 The splitting approximation over one step is therefore
 
 !bt
-\[ E = 1 - 4F\sin^p -q - (1-q)(1-4F\sin^2 p) = -q(2 - F\sin^2 p)) \]
+\[ E = 1 - 4F\sin^p -q - (1-q)(1-4F\sin^2 p) = -q(2 - F\sin^2 p))\tp \]
 !et