Exercises

Problem 1: Determine if equations are nonlinear or not

Classify each term in the following equations as linear or nonlinear. Assume that \(a\) and \(b\) are unknown numbers and that \(u\) and \(v\) are unknown functions. All other symbols are known quantities.

1. \(b^2 = 1\)
2. \(a+b = 1\), \(a-2b = 0\)
3. \(mu'' + \beta |u'|u' + cu = F(t)\)
4. \(u_t = {\alpha} u_{xx}\)
5. \(u_{tt} = c^2\nabla^2 u\)
6. \(u_t = \nabla\cdot({\alpha}(u)\nabla u) + f(x,y)\)
7. \(u_t + f(u)_x = 0\)
8. \(\boldsymbol{u}_t + \boldsymbol{u}\cdot\nabla \boldsymbol{u} = -\nabla p + r\nabla^2\boldsymbol{u}\), \(\nabla\cdot\boldsymbol{u} = 0\) (\(\boldsymbol{u}\) is a vector field)
9. \(u' = f(u,t)\)
10. \(\nabla^2 u = \lambda e^u\)

Problem 2: Linearize a nonlinear vibration ODE

Consider a nonlinear vibration problem

\[mu'' + bu'|u'| + s(u) = F(t),\]

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''\) is mass times acceleration, and \(bu'|u'|\) models water or air drag.

Approximate \(u''\) by a centered finite difference \(D_tD_t u\), and use a centered difference \(D_t u\) for \(u'\) as well. Observe then that \(s(u)\) does not contribute to making the resulting algebraic equation at a time level nonlinear. Use a geometric mean to linearize the quadratic nonlinearity arising from the term \(bu'|u'|\).

Exercise 3: Find the sparsity of the Jacobian

Consider a typical nonlinear Laplace term like \(\nabla\cdot{\alpha}(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 \({\alpha}\nabla^2 u\).

Hint. Set up the unknowns that enter the difference stencil and find the sparsity of the Jacobian that arise from the stencil.

Filename: nonlin_sparsity_Jacobian.pdf.

Exercise 4: 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: nonlin_Newton_linear.pdf.

Exercise 5: Differentiate a highly nonlinear term

The operator \(\nabla\cdot({\alpha}(u)\nabla u)\) with \({\alpha}(u) = ||\nabla u||^q\) appears in several physical problems, especially flow of Non-Newtonian fluids. In a Newton method one has to carry out the differentiation \(\partial{\alpha}(u)/\partial c_j\), for \(u=\sum_kc_k{\psi}_k\). Show that

\[{\partial\over\partial u_j} ||\nabla u||^q = q||\nabla u||^{q-2}\nabla u\cdot \nabla{\psi}_j{\thinspace .}\]

Filename: nonlin_differentiate.pdf.

Problem 6: Discretize a 1D problem with a nonlinear coefficient

We consider the problem

(1)\[ ((1 + u^2)u')' = 1,\quad x\in (0,1),\quad u(0)=u(1)=0{\thinspace .}\]

a) Discretize (1) by a centered finite difference method on a uniform mesh.

b) Discretize (1) by a finite element method with P1 of equal length. Use the Trapezoidal method to compute all integrals. Set up the resulting matrix system.

Filename: nonlin_1D_coeff_discretize.pdf.

Problem 7: Linearize a 1D problem with a nonlinear coefficient

We have a two-point boundary value problem

(2)\[ ((1 + u^2)u')' = 1,\quad x\in (0,1),\quad u(0)=u(1)=0{\thinspace .}\]

a) Construct a Picard iteration method for (2) without discretizing in space.

b) Apply Newton’s method to (2) without discretizing in space.

c) Discretize (2) by a centered finite difference scheme. Construct a Picard method for the resulting system of nonlinear algebraic equations.

d) Discretize (2) 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.pdf.

Problem 8: Finite differences for the 1D Bratu problem

We address the so-called Bratu problem

(3)\[ u'' + \lambda e^u=0,\quad x\in (0,1),\quad u(0)=u(1)=0,\]

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 (3) has an exact solution

\[u(x) = -2\ln\left(\frac{\cosh((x-\frac{1}{2})\theta/2)}{\cosh(\theta/4)}\right),\]

where \(\theta\) solves

\[\theta = \sqrt{2\lambda}\cosh(\theta/4){\thinspace .}\]

There are two solutions of (3) 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){\thinspace .}\]

A numerical value is \(\lambda_c = 3.513830719\).

a) Discretize (3) 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.

Filename: nonlin_1D_Bratu_fd.pdf.

Problem 9: Integrate functions of finite element expansions

We shall investigate integrals on the form

(4)\[ \int_0^L f(\sum_ku_k{\varphi}_k(x)){\varphi}_i(x){\, \mathrm{d}x},\]

where \({\varphi}_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 \(x_{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 \([x_{i-1},x_{i+1}]\). On this interval \({\varphi}_k(x)=0\) for \(k<i-1\) and \(k>i+1\), so only three basis functions will contribute:

\[\sum_k u_k{\varphi}_k(x) = u_{i-1}{\varphi}_{i-1}(x) + u_{i}{\varphi}_{i}(x) + u_{i+1}{\varphi}_{i+1}(x){\thinspace .}\]

The integral (4) now takes the simplified form

\[\int_{x_{i-1}}^{x_{i+1}} f(u_{i-1}{\varphi}_{i-1}(x) + u_{i}{\varphi}_{i}(x) + u_{i+1}{\varphi}_{i+1}(x)){\varphi}_i(x){\, \mathrm{d}x}{\thinspace .}\]

Split this integral in two integrals over cell L (left), \([x_{i-1},x_{i}]\), and cell R (right), \([x_{i},x_{i+1}]\). Over cell L, \(u\) simplifies to \(u_{i-1}{\varphi}_{i-1} + u_{i}{\varphi}_{i}\) (since \({\varphi}_{i+1}=0\) on this cell), and over cell R, \(u\) simplifies to \(u_{i}{\varphi}_{i} + u_{i+1}{\varphi}_{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 make LaTeX code and render it via http://latex.codecogs.com. 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.py.

Problem 10: Finite elements for the 1D Bratu problem

We address the same 1D Bratu problem as described in Problem 8: Finite differences for the 1D Bratu problem.

a) Discretize (Problem 10: 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.pdf.

Problem 11: Derive the Newton system 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.

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_Newton.pdf.

Problem 12: Derive algebraic equations for nonlinear 1D heat conduction

Consider a 1D heat conduction PDE

\[\varrho c(T) T_t = (k(T)T_x)_x,\]

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.

Use a uniform finite element mesh, P1 elements, and the group finite element method to derive the algebraic equations arising from the heat conduction PDE

a) Discretize the PDE by a finite difference method. Use either a Backward Euler or Crank-Nicolson scheme in time.

b) Derive the matrix and right-hand side of a Newton method applied to the discretized PDE.

Filename: nonlin_1D_heat_PDE.pdf.

Problem 13: 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'(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.

Bibliography

[Ref1]

C. T. Kelley. Iterative Methods for Linear and Nonlinear Equations, SIAM, 1995.

[Ref2]

M. Mortensen, H. P. Langtangen and G. N. Wells. A FEniCS-Based Programming Framework for Modeling Turbulent Flow by the Reynolds-Averaged Navier-Stokes Equations, Advances in Water Resources, 34(9), doi: 10.1016/j.advwatres.2011.02.013, 2011.