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{92} \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{93} \end{equation} $$ where \( p>0 \) is some real number.
a) Formulate a Forward Euler, Backward Euler, and a Crank-Nicolson scheme for (92).
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 (93).
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 models (93) and \eqref{nonlin:exer:logistic:gen:r2}. 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 (18)-(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{94} \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{95} \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{96} \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 (95) and (96) 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 (95) and (96) 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
(95) and
(96) 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{97} \end{equation} $$
a) Discretize (97) by a centered finite difference method on a uniform mesh.
b) Discretize (97) by a finite element method with P1 elements of equal length. Use the Trapezoidal method to compute all integrals. Set up the resulting matrix system in symbolic form such that the equations can be compared with those in a).
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{98} \end{equation} $$
a) Construct a Picard iteration method for (98) without discretizing in space.
b) Apply Newton's method to (98) without discretizing in space.
c) Discretize (98) by a centered finite difference scheme. Construct a Picard method for the resulting system of nonlinear algebraic equations.
d) Discretize (98) 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{99} \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 (99) 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 (99) 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 (99) 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 shall investigate integrals on the form
$$
\begin{equation}
\int_0^L f(\sum_ku_k\basphi_k(x))\basphi_i(x)\dx,
\tag{100}
\end{equation}
$$
where \( \basphi_i(x) \) are P1 finite element basis functions and \( u_k \)
are unknown coefficients, more precisely the values of the unknown
function \( u \) at nodes \( \xno{k} \). We introduce a node numbering that
goes from left to right and also that all cells have
the same length \( h \). Given \( i \), the integral
only gets contributions from \( [\xno{i-1},\xno{i+1}] \). On this
interval \( \basphi_k(x)=0 \) for \( k < i-1 \) and \( k>i+1 \), so only three
basis functions will contribute:
$$
\sum_k u_k\basphi_k(x) = u_{i-1}\basphi_{i-1}(x) +
u_{i}\basphi_{i}(x) + u_{i+1}\basphi_{i+1}(x)\tp
$$
The integral (100) now takes the
simplified form
$$
\int_{\xno{i-1}}^{\xno{i+1}}
f(u_{i-1}\basphi_{i-1}(x) +
u_{i}\basphi_{i}(x) + u_{i+1}\basphi_{i+1}(x))\basphi_i(x)\dx\tp
$$
Split this integral in two integrals over cell L (left),
\( [\xno{i-1},\xno{i}] \), and cell R (right), \( [\xno{i},\xno{i+1}] \). Over
cell L, \( u \) simplifies to \( u_{i-1}\basphi_{i-1} + u_{i}\basphi_{i} \)
(since \( \basphi_{i+1}=0 \) on this cell), and over cell R, \( u \)
simplifies to \( u_{i}\basphi_{i} + u_{i+1}\basphi_{i+1} \). Make a
sympy program that can compute the integral and write it out as a
difference equation. Give the \( f(u) \) formula on the command line.
Try out \( f(u)=u^2, \sin u, \exp u \).
Introduce symbols u_i, u_im1, and u_ip1 for \( u_i \), \( u_{i-1} \),
and \( u_{i+1} \), respectively, and similar symbols for \( x_i \), \( x_{i-1} \),
and \( x_{i+1} \). Find formulas for the basis functions on each of the
two cells, make expressions for \( u \) on the two cells, integrate over
each cell, expand the answer and simplify. You can ask sympy for
LaTeX code and render it either by creating a LaTeX document and
compiling it to a PDF document or by using
http://latex.codecogs.com to display LaTeX formulas in a web
page. Here are some appropriate Python statements
for the latter purpose:
from sympy import *
...
# expr_i holdes the integral as a sympy expression
latex_code = latex(expr_i, mode='plain')
# Replace u_im1 sympy symbol name by latex symbol u_{i-1}
latex_code = latex_code.replace('im1', '{i-1}')
# Replace u_ip1 sympy symbol name by latex symbol u_{i+1}
latex_code = latex_code.replace('ip1', '{i+1}')
# Escape (quote) latex_code so it can be sent as HTML text
import cgi
html_code = cgi.escape(latex_code)
# Make a file with HTML code for displaying the LaTeX formula
f = open('tmp.html', 'w')
# Include an image that can be clicked on to yield a new
# page with an interactive editor and display area where the
# formula can be further edited
text = """
<a href="http://www.codecogs.com/eqnedit.php?latex=%(html_code)s"
target="_blank">
<img src="http://latex.codecogs.com/gif.latex?%(html_code)s"
title="%(latex_code)s"/>
</a>
""" % vars()
f.write(text)
f.close()
The formula is displayed by loading tmp.html into a web browser.
Filename: fu_fem_int.
We address the same 1D Bratu problem as described in Problem 10: Finite differences for the 1D Bratu problem.
a) Discretize (Problem 12: Finite elements for the 1D Bratu problem) by a finite element method using a uniform mesh with P1 elements. Use a group finite element method for the \( e^u \) term.
b) Set up the nonlinear equations \( F_i(u_0,u_1,\ldots,u_{N_x})=0 \) from a). Calculate the associated Jacobian.
Filename: nonlin_1D_Bratu_fe.
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 symbol \( u \) has several meanings, depending on the context, as briefly mentioned in the section Finite element discretization. Go through the derivation of the Picard iteration method in that section and use different symbols for all the different approximations of \( u \):
nonlin_heat_FE_usymbols.
We study the multi-dimensional heat conduction PDE $$ \varrho c(T) T_t = \nabla\cdot (k(T)\nabla T)$$ in a spatial domain \( \Omega \), with a nonlinear Robin boundary condition $$ -k(T)\frac{\partial T}{\partial n} = h(T)(T-T_s(t)),$$ at the boundary \( \partial\Omega \). The primary unknown is the temperature \( T \), \( \varrho \) is the density of the solid material, \( c(T) \) is the heat capacity, \( k(T) \) is the heat conduction, \( h(T) \) is a heat transfer coefficient, and \( T_s(T) \) is a possibly time-dependent temperature of the surroundings.
a) Use a Backward Euler or Crank-Nicolson time discretization and derive the variational form for the spatial problem to be solved at each time level.
b) Define a Picard iteration method from the variational form at a time level.
c) Derive expressions for the matrix and the right-hand side of the equation system that arises from applying Newton's method to the variational form at a time level.
d) Apply the Backward Euler or Crank-Nicolson scheme in time first. Derive a Newton method at the PDE level. Make a variational form of the resulting PDE at a time level.
Filename: nonlin_heat_FE.
We consider the same problem as in Exercise 15: Derive Picard and Newton systems from a variational form, but restricted to one space dimension: \( \Omega = [0,L] \).
Simplify the boundary condition to \( T_x=0 \) (i.e., \( h(T)=0 \)).
Use a uniform finite element mesh of P1 elements, the group
finite element method, and the Trapezoidal
rule for integration at the nodes to derive symbolic expressions for
the algebraic equations arising from this diffusion problem.
Filename: nonlin_1D_heat_FE.
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.