.. !split .. _nonlin:exer: Exercises ========= .. --- begin exercise --- .. _nonlin:exer:lin:vs:nonlin: Problem 1: Determine if equations are nonlinear or not ------------------------------------------------------ Classify each term in the following equations as linear or nonlinear. Assume that :math:`a` and :math:`b` are unknown numbers and that :math:`u` and :math:`v` are unknown functions. All other symbols are known quantities. 1. :math:`b^2 = 1` 2. :math:`a+b = 1`, :math:`a-2b = 0` 3. :math:`mu'' + \beta |u'|u' + cu = F(t)` 4. :math:`u_t = {\alpha} u_{xx}` 5. :math:`u_{tt} = c^2\nabla^2 u` 6. :math:`u_t = \nabla\cdot({\alpha}(u)\nabla u) + f(x,y)` 7. :math:`u_t + f(u)_x = 0` 8. :math:`\boldsymbol{u}_t + \boldsymbol{u}\cdot\nabla \boldsymbol{u} = -\nabla p + r\nabla^2\boldsymbol{u}`, :math:`\nabla\cdot\boldsymbol{u} = 0` (:math:`\boldsymbol{u}` is a vector field) 9. :math:`u' = f(u,t)` 10. :math:`\nabla^2 u = \lambda e^u` .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:vib:geometric:mean: Problem 2: Linearize a nonlinear vibration ODE ---------------------------------------------- Consider a nonlinear vibration problem .. math:: mu'' + bu'|u'| + s(u) = F(t), where :math:`m>0` is a constant, :math:`b\geq 0` is a constant, :math:`s(u)` a possibly nonlinear function of :math:`u`, and :math:`F(t)` is a prescribed function. Such models arise from Newton's second law of motion in mechanical vibration problems where :math:`s(u)` is a spring or restoring force, :math:`mu''` is mass times acceleration, and :math:`bu'|u'|` models water or air drag. Approximate :math:`u''` by a centered finite difference :math:`D_tD_t u`, and use a centered difference :math:`D_t u` for :math:`u'` as well. Observe then that :math:`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 :math:`bu'|u'|`. .. 2DO: b) Newmark scheme .. derive it logically and connect it to the centered diff scheme .. ma + bv|v| + s(u) = F(t), v'=a, u'=v (staggered is natural, .. v at n+1/2 and a and u at n). Should be in vib first .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:sparsity:Jacobian: Exercise 3: Find the sparsity of the Jacobian --------------------------------------------- Consider a typical nonlinear Laplace term like :math:`\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 :math:`{\alpha}\nabla^2 u`. .. --- begin hint in exercise --- **Hint.** Set up the unknowns that enter the difference stencil and find the sparsity of the Jacobian that arise from the stencil. .. --- end hint in exercise --- Filename: ``nonlin_sparsity_Jacobian.pdf``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:Newton:linear: Exercise 4: Newton's method for linear problems ----------------------------------------------- Suppose we have a linear system :math:`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``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:grad:pow:term: Exercise 5: Differentiate a highly nonlinear term ------------------------------------------------- The operator :math:`\nabla\cdot({\alpha}(u)\nabla u)` with :math:`{\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 :math:`\partial{\alpha}(u)/\partial c_j`, for :math:`u=\sum_kc_k{\psi}_k`. Show that .. math:: {\partial\over\partial u_j} ||\nabla u||^q = q||\nabla u||^{q-2}\nabla u\cdot \nabla{\psi}_j{\thinspace .} Filename: ``nonlin_differentiate.pdf``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:1pu2:fem: Problem 6: Discretize a 1D problem with a nonlinear coefficient --------------------------------------------------------------- We consider the problem .. _Eq:nonlin:exer:1D:1pu2:fem:pde: .. math:: :label: nonlin:exer:1D:1pu2:fem:pde ((1 + u^2)u')' = 1,\quad x\in (0,1),\quad u(0)=u(1)=0{\thinspace .} **a)** Discretize :eq:`nonlin:exer:1D:1pu2:fem:pde` by a centered finite difference method on a uniform mesh. **b)** Discretize :eq:`nonlin:exer:1D:1pu2:fem:pde` 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``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:1pu2:PicardNewton: Problem 7: Linearize a 1D problem with a nonlinear coefficient -------------------------------------------------------------- We have a two-point boundary value problem .. _Eq:nonlin:exer:1D:1pu2:PicardNewton:pde: .. math:: :label: nonlin:exer:1D:1pu2:PicardNewton:pde ((1 + u^2)u')' = 1,\quad x\in (0,1),\quad u(0)=u(1)=0{\thinspace .} **a)** Construct a Picard iteration method for :eq:`nonlin:exer:1D:1pu2:PicardNewton:pde` without discretizing in space. **b)** Apply Newton's method to :eq:`nonlin:exer:1D:1pu2:PicardNewton:pde` without discretizing in space. **c)** Discretize :eq:`nonlin:exer:1D:1pu2:PicardNewton:pde` by a centered finite difference scheme. Construct a Picard method for the resulting system of nonlinear algebraic equations. **d)** Discretize :eq:`nonlin:exer:1D:1pu2:PicardNewton:pde` 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``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:fu:discretize:fd: Problem 8: Finite differences for the 1D Bratu problem ------------------------------------------------------ We address the so-called Bratu problem .. _Eq:nonlin:exer:1D:fu:discretize:fd:pde: .. math:: :label: nonlin:exer:1D:fu:discretize:fd:pde u'' + \lambda e^u=0,\quad x\in (0,1),\quad u(0)=u(1)=0, where :math:`\lambda` is a given parameter and :math:`u` is a function of :math:`x`. This is a widely used model problem for studying numerical methods for nonlinear differential equations. The problem :eq:`nonlin:exer:1D:fu:discretize:fd:pde` has an exact solution .. math:: u(x) = -2\ln\left(\frac{\cosh((x-\frac{1}{2})\theta/2)}{\cosh(\theta/4)}\right), where :math:`\theta` solves .. math:: \theta = \sqrt{2\lambda}\cosh(\theta/4){\thinspace .} There are two solutions of :eq:`nonlin:exer:1D:fu:discretize:fd:pde` for :math:`0<\lambda <\lambda_c` and no solution for :math:`\lambda >\lambda_c`. For :math:`\lambda = \lambda_c` there is one unique solution. The critical value :math:`\lambda_c` solves .. math:: 1 = \sqrt{2\lambda_c}\frac{1}{4}\sinh(\theta(\lambda_c)/4){\thinspace .} A numerical value is :math:`\lambda_c = 3.513830719`. **a)** Discretize :eq:`nonlin:exer:1D:fu:discretize:fd:pde` by a centered finite difference method. **b)** Set up the nonlinear equations :math:`F_i(u_0,u_1,\ldots,u_{N_x})=0` from a). Calculate the associated Jacobian. Filename: ``nonlin_1D_Bratu_fd.pdf``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:fu:fem:int: Problem 9: Integrate functions of finite element expansions ----------------------------------------------------------- .. index:: latex.codecogs.com web site .. index:: online rendering of LaTeX formulas We shall investigate integrals on the form .. _Eq:nonlin:exer:fu:fem:int:global: .. math:: :label: nonlin:exer:fu:fem:int:global \int_0^L f(\sum_ku_k{\varphi}_k(x)){\varphi}_i(x){\, \mathrm{d}x}, where :math:`{\varphi}_i(x)` are P1 finite element basis functions and :math:`u_k` are unknown coefficients, more precisely the values of the unknown function :math:`u` at nodes :math:`x_{k}`. We introduce a node numbering that goes from left to right and also that all cells have the same length :math:`h`. Given :math:`i`, the integral only gets contributions from :math:`[x_{i-1},x_{i+1}]`. On this interval :math:`{\varphi}_k(x)=0` for :math:`ki+1`, so only three basis functions will contribute: .. math:: \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 :eq:`nonlin:exer:fu:fem:int:global` now takes the simplified form .. math:: \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), :math:`[x_{i-1},x_{i}]`, and cell R (right), :math:`[x_{i},x_{i+1}]`. Over cell L, :math:`u` simplifies to :math:`u_{i-1}{\varphi}_{i-1} + u_{i}{\varphi}_{i}` (since :math:`{\varphi}_{i+1}=0` on this cell), and over cell R, :math:`u` simplifies to :math:`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 :math:`f(u)` formula on the command line. Try out :math:`f(u)=u^2, \sin u, \exp u`. .. --- begin hint in exercise --- **Hint.** Introduce symbols ``u_i``, ``u_im1``, and ``u_ip1`` for :math:`u_i`, :math:`u_{i-1}`, and :math:`u_{i+1}`, respectively, and similar symbols for :math:`x_i`, :math:`x_{i-1}`, and :math:`x_{i+1}`. Find formulas for the basis functions on each of the two cells, make expressions for :math:`u` on the two cells, integrate over each cell, expand the answer and simplify. You can make LaTeX code and render it via ``_. Here are some appropriate Python statements for the latter purpose: .. code-block:: python 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 = """ """ % vars() f.write(text) f.close() The formula is displayed by loading ``tmp.html`` into a web browser. .. --- end hint in exercise --- Filename: ``fu_fem_int.py``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:fu:discretize:fe: Problem 10: Finite elements for the 1D Bratu problem ---------------------------------------------------- We address the same 1D Bratu problem as described in :ref:`nonlin:exer:1D:fu:discretize:fd`. **a)** Discretize (:ref:`nonlin:exer:1D:fu:discretize:fe`) by a finite element method using a uniform mesh with P1 elements. Use a group finite element method for the :math:`e^u` term. **b)** Set up the nonlinear equations :math:`F_i(u_0,u_1,\ldots,u_{N_x})=0` from a). Calculate the associated Jacobian. Filename: ``nonlin_1D_Bratu_fe.pdf``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:dD:heat:nonlinear:c:a: Problem 11: Derive the Newton system from a variational form ------------------------------------------------------------ We study the multi-dimensional heat conduction PDE .. math:: \varrho c(T) T_t = \nabla\cdot (k(T)\nabla T) in a spatial domain :math:`\Omega`, with a nonlinear Robin boundary condition .. math:: -k(T)\frac{\partial T}{\partial n} = h(T)(T-T_s(t)), at the boundary :math:`\partial\Omega`. The primary unknown is the temperature :math:`T`, :math:`\varrho` is the density of the solid material, :math:`c(T)` is the heat capacity, :math:`k(T)` is the heat conduction, :math:`h(T)` is a heat transfer coefficient, and :math:`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``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:heat:nonlinear:c:a: Problem 12: Derive algebraic equations for nonlinear 1D heat conduction ----------------------------------------------------------------------- Consider a 1D heat conduction PDE .. math:: \varrho c(T) T_t = (k(T)T_x)_x, where :math:`\varrho` is the density of the solid material, :math:`c(T)` is the heat capacity, :math:`T` is the temperature, and :math:`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``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:continuation:1DnNflow: Problem 13: Investigate a 1D problem with a continuation method --------------------------------------------------------------- .. index:: continuation method Flow of a pseudo-plastic power-law fluid between two flat plates can be modeled by .. math:: \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 :math:`\beta>0` and :math:`\mu_0>0` are constants. A target value of :math:`n` may be :math:`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 :math:`u` on a mesh. Verify that the solver gives an exact solution for :math:`n=1` on a uniform mesh regardless of the cell size. **c)** Given a sequence of decreasing :math:`n` values, solve the problem for each :math:`n` using the solution for the previous :math:`n` as initial guess for the Picard iteration. This is called a continuation method. Experiment with :math:`n=(1,0.6,0.2)` and :math:`n=(1,0.9,0.8,\ldots,0.2)` and make a table of the number of Picard iterations versus :math:`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. .. --- end exercise --- 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.