Classify each term in the following equations as linear or nonlinear. Assume that \( u \), \( \u \), and \( p \) are unknown functions and that all other symbols are known quantities.
nonlinear_vs_linear.
The logistic model for population growth is derived by assuming a nonlinear growth rate, $$ \begin{equation} u^{\prime} = a(u)u,\quad u(0)=I, \tag{5.82} \end{equation} $$ and the logistic model arises from the simplest possible choice of \( a(u) \): \( r(u)=\varrho(1 - u/M) \), where \( M \) is the maximum value of \( u \) that the environment can sustain, and \( \varrho \) is the growth under unlimited access to resources (as in the beginning when \( u \) is small). The idea is that \( a(u)\sim\varrho \) when \( u \) is small and that \( a(t)\rightarrow 0 \) as \( u\rightarrow M \).
An \( a(u) \) that generalizes the linear choice is the polynomial form $$ \begin{equation} a(u) = \varrho(1-u/M)^p, \tag{5.83} \end{equation} $$ where \( p>0 \) is some real number.
a) Formulate a Forward Euler, Backward Euler, and a Crank-Nicolson scheme for (5.82).
Use a geometric mean approximation in the Crank-Nicolson scheme: \( [a(u)u]^{n+1/2}\approx a(u^n)u^{n+1} \).
b) Formulate Picard and Newton iteration for the Backward Euler scheme in a).
c) Implement the numerical solution methods from a) and b). Use logistic.py to compare the case \( p=1 \) and the choice (5.83).
d) Implement unit tests that check the asymptotic limit of the solutions: \( u\rightarrow M \) as \( t\rightarrow\infty \).
You need to experiment to find what "infinite time" is (increases substantially with \( p \)) and what the appropriate tolerance is for testing the asymptotic limit.
e) Perform experiments with Newton and Picard iteration for the model (5.83). See how sensitive the number of iterations is to \( \Delta t \) and \( p \).
Filename: logistic_p.
The program Newton_demo.py illustrates graphically each step in Newton's method and is run like
Terminal> python Newton_demo.py f dfdx x0 xmin xmax
Use this program to investigate potential problems with Newton's method when solving \( e^{-0.5x^2}\cos (\pi x)=0 \). Try a starting point \( x_0=0.8 \) and \( x_0=0.85 \) and watch the different behavior. Just run
Terminal> python Newton_demo.py '0.2 + exp(-0.5*x**2)*cos(pi*x)' \
'-x*exp(-x**2)*cos(pi*x) - pi*exp(-x**2)*sin(pi*x)' \
0.85 -3 3
and repeat with 0.85 replaced by 0.8.
Write up the system (5.18)-(5.19) in the form \( F(u)=0 \), \( F=(F_0,F_1) \), \( u=(u_0,u_1) \), and compute the Jacobian \( J_{i,j}=\partial F_i/\partial u_j \).
Consider a nonlinear vibration problem $$ \begin{equation} mu^{\prime\prime} + bu^{\prime}|u^{\prime}| + s(u) = F(t), \tag{5.84} \end{equation} $$ where \( m>0 \) is a constant, \( b\geq 0 \) is a constant, \( s(u) \) a possibly nonlinear function of \( u \), and \( F(t) \) is a prescribed function. Such models arise from Newton's second law of motion in mechanical vibration problems where \( s(u) \) is a spring or restoring force, \( mu^{\prime\prime} \) is mass times acceleration, and \( bu^{\prime}|u^{\prime}| \) models water or air drag.
a) Rewrite the equation for \( u \) as a system of two first-order ODEs, and discretize this system by a Crank-Nicolson (centered difference) method. With \( v=u^\prime \), we get a nonlinear term \( v^{n+\frac{1}{2}}|v^{n+\frac{1}{2}}| \). Use a geometric average for \( v^{n+\frac{1}{2}} \).
b) Formulate a Picard iteration method to solve the system of nonlinear algebraic equations.
c) Explain how to apply Newton's method to solve the nonlinear equations at each time level. Derive expressions for the Jacobian and the right-hand side in each Newton iteration.
Filename: nonlin_vib.
In the section Crank-Nicolson discretization we introduce alternative arithmetic means of a product. Say the product is \( P(t)Q(t) \) evaluated at \( t=t_{n+\frac{1}{2}} \). The exact value is $$ [PQ]^{n+\frac{1}{2}} = P^{n+\frac{1}{2}}Q^{n+\frac{1}{2}} $$ There are two obvious candidates for evaluating \( [PQ]^{n+\frac{1}{2}} \) as a mean of values of \( P \) and \( Q \) at \( t_n \) and \( t_{n+1} \). Either we can take the arithmetic mean of each factor \( P \) and \( Q \), $$ \begin{equation} [PQ]^{n+\frac{1}{2}} \approx \frac{1}{2}(P^n + P^{n+1})\frac{1}{2}(Q^n + Q^{n+1}), \tag{5.85} \end{equation} $$ or we can take the arithmetic mean of the product \( PQ \): $$ \begin{equation} [PQ]^{n+\frac{1}{2}} \approx \frac{1}{2}(P^nQ^n + P^{n+1}Q^{n+1})\tp \tag{5.86} \end{equation} $$
The arithmetic average of \( P(t_{n+\frac{1}{2}}) \) is \( \Oof{\Delta t^2} \): $$ P(t_{n+\frac{1}{2}}) = \frac{1}{2}(P^n + P^{n+1}) +\Oof{\Delta t^2}\tp$$ A fundamental question is whether (5.85) and (5.86) have different orders of accuracy in \( \Delta t = t_{n+1}-t_n \). To investigate this question, expand quantities at \( t_{n+1} \) and \( t_n \) in Taylor series around \( t_{n+\frac{1}{2}} \), and subtract the true value \( [PQ]^{n+\frac{1}{2}} \) from the approximations (5.85) and (5.86) to see what the order of the error terms are.
You may explore sympy for carrying out the tedious calculations.
A general Taylor series expansion of \( P(t+\frac{1}{2}\Delta t) \) around \( t \)
involving just a general function \( P(t) \) can be
created as follows:
>>> from sympy import *
>>> t, dt = symbols('t dt')
>>> P = symbols('P', cls=Function)
>>> P(t).series(t, 0, 4)
P(0) + t*Subs(Derivative(P(_x), _x), (_x,), (0,)) +
t**2*Subs(Derivative(P(_x), _x, _x), (_x,), (0,))/2 +
t**3*Subs(Derivative(P(_x), _x, _x, _x), (_x,), (0,))/6 + O(t**4)
>>> P_p = P(t).series(t, 0, 4).subs(t, dt/2)
>>> P_p
P(0) + dt*Subs(Derivative(P(_x), _x), (_x,), (0,))/2 +
dt**2*Subs(Derivative(P(_x), _x, _x), (_x,), (0,))/8 +
dt**3*Subs(Derivative(P(_x), _x, _x, _x), (_x,), (0,))/48 + O(dt**4)
The error of the arithmetic mean, \( \frac{1}{2}(P(-\frac{1}{2}\Delta t) + P(-\frac{1}{2}\Delta t)) \) for \( t=0 \) is then
>>> P_m = P(t).series(t, 0, 4).subs(t, -dt/2)
>>> mean = Rational(1,2)*(P_m + P_p)
>>> error = simplify(expand(mean) - P(0))
>>> error
dt**2*Subs(Derivative(P(_x), _x, _x), (_x,), (0,))/8 + O(dt**4)
Use these examples to investigate the error of
(5.85) and
(5.86) for \( n=0 \). (Choosing \( n=0 \)
is necessary for not making the expressions too complicated for sympy,
but there is of course no lack of generality by using \( n=0 \) rather
than an arbitrary \( n \) - the main point is the product and addition
of Taylor series.)
Filename: product_arith_mean.
Suppose we have a linear system \( F(u) = Au- b=0 \). Apply Newton's method
to this system, and show that the method converges in one iteration.
Filename: Newton_linear.
We consider the problem $$ \begin{equation} ((1 + u^2)u^{\prime})^{\prime} = 1,\quad x\in (0,1),\quad u(0)=u(1)=0\tp \tag{5.87} \end{equation} $$
Discretize (5.87) by a centered
finite difference method on a uniform mesh.
Filename: nonlin_1D_coeff_discretize.
We have a two-point boundary value problem $$ \begin{equation} ((1 + u^2)u^{\prime})^{\prime} = 1,\quad x\in (0,1),\quad u(0)=u(1)=0\tp \tag{5.88} \end{equation} $$
a) Construct a Picard iteration method for (5.88) without discretizing in space.
b) Apply Newton's method to (5.88) without discretizing in space.
c) Discretize (5.88) by a centered finite difference scheme. Construct a Picard method for the resulting system of nonlinear algebraic equations.
d) Discretize (5.88) by a centered finite difference scheme. Define the system of nonlinear algebraic equations, calculate the Jacobian, and set up Newton's method for solving the system.
Filename: nonlin_1D_coeff_linearize.
We address the so-called Bratu problem $$ \begin{equation} u^{\prime\prime} + \lambda e^u=0,\quad x\in (0,1),\quad u(0)=u(1)=0, \tag{5.89} \end{equation} $$ where \( \lambda \) is a given parameter and \( u \) is a function of \( x \). This is a widely used model problem for studying numerical methods for nonlinear differential equations. The problem (5.89) has an exact solution $$ \uex(x) = -2\ln\left(\frac{\cosh((x-\half)\theta/2)}{\cosh(\theta/4)}\right),$$ where \( \theta \) solves $$ \theta = \sqrt{2\lambda}\cosh(\theta/4)\tp$$ There are two solutions of (5.89) for \( 0 < \lambda < \lambda_c \) and no solution for \( \lambda >\lambda_c \). For \( \lambda = \lambda_c \) there is one unique solution. The critical value \( \lambda_c \) solves $$ 1 = \sqrt{2\lambda_c}\frac{1}{4}\sinh(\theta(\lambda_c)/4)\tp$$ A numerical value is \( \lambda_c = 3.513830719 \).
a) Discretize (5.89) by a centered finite difference method.
b) Set up the nonlinear equations \( F_i(u_0,u_1,\ldots,u_{N_x})=0 \) from a). Calculate the associated Jacobian.
c) Implement a solver that can compute \( u(x) \) using Newton's method. Plot the error as a function of \( x \) in each iteration.
d) Investigate whether Newton's method gives second-order convergence by computing \( || \uex - u||/||\uex - u^{-}||^2 \) in each iteration, where \( u \) is solution in the current iteration and \( u^{-} \) is the solution in the previous iteration.
Filename: nonlin_1D_Bratu_fd.
We address the 1D heat conduction PDE $$ \varrho c(T) T_t = (k(T)T_x)_x,$$ for \( x\in [0,L] \), where \( \varrho \) is the density of the solid material, \( c(T) \) is the heat capacity, \( T \) is the temperature, and \( k(T) \) is the heat conduction coefficient. \( T(x,0)=I(x) \), and ends are subject to a cooling law: $$ k(T)T_x|_{x=0} = h(T)(T-T_s),\quad -k(T)T_x|_{x=L}=h(T)(T-T_s),$$ where \( h(T) \) is a heat transfer coefficient and \( T_s \) is the given surrounding temperature.
a) Discretize this PDE in time using either a Backward Euler or Crank-Nicolson scheme.
b) Formulate a Picard iteration method for the time-discrete problem (i.e., an iteration method before discretizing in space).
c) Formulate a Newton method for the time-discrete problem in b).
d) Discretize the PDE by a finite difference method in space. Derive the matrix and right-hand side of a Picard iteration method applied to the space-time discretized PDE.
e) Derive the matrix and right-hand side of a Newton method applied to the discretized PDE in d).
Filename: nonlin_1D_heat_FD.
The operator \( \nabla\cdot(\dfc(u)\nabla u) \) with \( \dfc(u) = |\nabla u|^q \) appears in several physical problems, especially flow of Non-Newtonian fluids. The expression \( |\nabla u| \) is defined as the Euclidean norm of a vector: \( |\nabla u|^2 = \nabla u \cdot \nabla u \). In a Newton method one has to carry out the differentiation \( \partial\dfc(u)/\partial c_j \), for \( u=\sum_kc_k\baspsi_k \). Show that $$ {\partial\over\partial u_j} |\nabla u|^q = q|\nabla u|^{q-2}\nabla u\cdot \nabla\baspsi_j\tp $$
Filename: nonlin_differentiate.
Redo the section Finite difference discretization when a Crank-Nicolson scheme is used to discretize the equations in time and the problem is formulated for three spatial dimensions.
Express the Jacobian as \( J_{i,j,k,r,s,t} = \partial F_{i,j,k}/\partial u_{r,s,t} \) and observe, as in the 2D case, that \( J_{i,j,k,r,s,t} \) is very sparse: \( J_{i,j,k,r,s,t}\neq 0 \) only for \( r=i\pm i \), \( s=j\pm 1 \), and \( t=k\pm 1 \) as well as \( r=i \), \( s=j \), and \( t=k \).
Filename: nonlin_heat_FD_CN_2D.
Consider a typical nonlinear Laplace term like \( \nabla\cdot\dfc(u)\nabla u \) discretized by centered finite differences. Explain why the Jacobian corresponding to this term has the same sparsity pattern as the matrix associated with the corresponding linear term \( \dfc\nabla^2 u \).
Set up the unknowns that enter the difference equation at a point \( (i,j) \) in 2D or \( (i,j,k) \) in 3D, and identify the nonzero entries of the Jacobian that can arise from such a type of difference equation.
Filename: nonlin_sparsity_Jacobian.
Flow of a pseudo-plastic power-law fluid between two flat plates can be modeled by $$ \frac{d}{dx}\left(\mu_0\left\vert\frac{du}{dx}\right\vert^{n-1} \frac{du}{dx}\right) = -\beta,\quad u^{\prime}(0)=0,\ u(H) = 0,$$ where \( \beta>0 \) and \( \mu_0>0 \) are constants. A target value of \( n \) may be \( n=0.2 \).
a) Formulate a Picard iteration method directly for the differential equation problem.
b) Perform a finite difference discretization of the problem in each Picard iteration. Implement a solver that can compute \( u \) on a mesh. Verify that the solver gives an exact solution for \( n=1 \) on a uniform mesh regardless of the cell size.
c) Given a sequence of decreasing \( n \) values, solve the problem for each \( n \) using the solution for the previous \( n \) as initial guess for the Picard iteration. This is called a continuation method. Experiment with \( n=(1,0.6,0.2) \) and \( n=(1,0.9,0.8,\ldots,0.2) \) and make a table of the number of Picard iterations versus \( n \).
d) Derive a Newton method at the differential equation level and discretize the resulting linear equations in each Newton iteration with the finite difference method.
e) Investigate if Newton's method has better convergence properties than Picard iteration, both in combination with a continuation method.