Solving nonlinear ODE and PDE problems

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.

Exercise 2: Derive and investigate a generalized logistic model

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).

Hint. Use a geometric mean approximation in the Crank-Nicolson scheme: $ [a(u)u]^{n+1/2}\approx a(u^n)u^{n+1} $.

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 $.

Hint. 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 $.

Problem 3: Experience the behavior of Newton's method

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

Problem 4: Compute the Jacobian of a $ 2\times 2 $ system

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 $.

Problem 5: Solve nonlinear equations arising from a vibration ODE

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.

Exercise 6: Find the truncation error of arithmetic mean of products

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.

Hint. 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.)

Problem 7: Newton's method for linear problems

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.

Exercise 8: Discretize a 1D problem with a nonlinear coefficient

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} $$

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).

Exercise 9: Linearize a 1D problem with a nonlinear coefficient

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} $$

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.

Problem 10: Finite differences for the 1D Bratu problem

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 $.

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.

Problem 11: Integrate functions of finite element expansions

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 $.

Hint. 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()

Problem 12: Finite elements for the 1D Bratu problem

b) Set up the nonlinear equations $ F_i(u_0,u_1,\ldots,u_{N_x})=0 $ from a). Calculate the associated Jacobian.

Exercise 13: Discretize a nonlinear 1D heat conduction PDE by finite differences

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).

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).

Exercise 14: Use different symbols for different approximations of the solution

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 $:

Exercise 15: Derive Picard and Newton systems from a variational form

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.

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.

Exercise 16: Derive algebraic equations for nonlinear 1D heat conduction

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.

Exercise 17: Differentiate a highly nonlinear term

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 $$

Exercise 18: Crank-Nicolson for a nonlinear 3D diffusion equation

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.

Hint. 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 $.

Exercise 19: Find the sparsity of the Jacobian

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 $.

Hint. 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.

Problem 20: Investigate a 1D problem with a continuation method

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.

Exercises

Problem 1: Determine if equations are nonlinear or not