.. !split .. _nonlin:exer: Exercises (9) ====================== .. --- begin exercise --- .. _nonlin:exer:lin:vs:nonlin: Problem 5.1: Determine if equations are nonlinear or not -------------------------------------------------------- Classify each term in the following equations as linear or nonlinear. Assume that :math:`u`, :math:`\boldsymbol{u}`, and :math:`p` are unknown functions and that all other symbols are known quantities. 1. :math:`mu^{\prime\prime} + \beta |u^{\prime}|u^{\prime} + cu = F(t)` 2. :math:`u_t = {\alpha} u_{xx}` 3. :math:`u_{tt} = c^2\nabla^2 u` 4. :math:`u_t = \nabla\cdot({\alpha}(u)\nabla u) + f(x,y)` 5. :math:`u_t + f(u)_x = 0` 6. :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) 7. :math:`u^{\prime} = f(u,t)` 8. :math:`\nabla^2 u = \lambda e^u` .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``nonlinear_vs_linear``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:logistic:gen: Problem 5.2: Derive and investigate a generalized logistic model ---------------------------------------------------------------- The logistic model for population growth is derived by assuming a nonlinear growth rate, .. _Eq:nonlin:exer:logistic:gen:eq: .. math:: \tag{601} u^{\prime} = a(u)u,\quad u(0)=I, and the logistic model arises from the simplest possible choice of :math:`a(u)`: :math:`r(u)=\varrho(1 - u/M)`, where :math:`M` is the maximum value of :math:`u` that the environment can sustain, and :math:`\varrho` is the growth under unlimited access to resources (as in the beginning when :math:`u` is small). The idea is that :math:`a(u)\sim\varrho` when :math:`u` is small and that :math:`a(t)\rightarrow 0` as :math:`u\rightarrow M`. An :math:`a(u)` that generalizes the linear choice is the polynomial form .. _Eq:nonlin:exer:logistic:gen:r1: .. math:: \tag{602} a(u) = \varrho(1-u/M)^p, where :math:`p>0` is some real number. **a)** Formulate a Forward Euler, Backward Euler, and a Crank-Nicolson scheme for :ref:`(601) `. .. --- begin hint in exercise --- **Hint.** Use a geometric mean approximation in the Crank-Nicolson scheme: :math:`[a(u)u]^{n+1/2}\approx a(u^n)u^{n+1}`. .. --- end hint in exercise --- .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) **b)** Formulate Picard and Newton iteration for the Backward Euler scheme in a). .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) **c)** Implement the numerical solution methods from a) and b). Use `logistic.py `__ to compare the case :math:`p=1` and the choice :ref:`(602) `. .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) **d)** Implement unit tests that check the asymptotic limit of the solutions: :math:`u\rightarrow M` as :math:`t\rightarrow\infty`. .. --- begin hint in exercise --- **Hint.** You need to experiment to find what "infinite time" is (increases substantially with :math:`p`) and what the appropriate tolerance is for testing the asymptotic limit. .. --- end hint in exercise --- .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) **e)** Perform experiments with Newton and Picard iteration for the model :ref:`(602) `. See how sensitive the number of iterations is to :math:`\Delta t` and :math:`p`. .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) .. ===== Exercise: Derive a relaxation formula ===== Filename: ``logistic_p``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:Newton:problems1: Problem 5.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 .. code-block:: text Terminal> python Newton_demo.py f dfdx x0 xmin xmax Use this program to investigate potential problems with Newton's method when solving :math:`e^{-0.5x^2}\cos (\pi x)=0`. Try a starting point :math:`x_0=0.8` and :math:`x_0=0.85` and watch the different behavior. Just run .. code-block:: text 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. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:vib:Jacobian: Exercise 5.4: Compute the Jacobian of a :math:`2\times 2` system ---------------------------------------------------------------- Write up the system :ref:`(537) `-:ref:`(538) ` in the form :math:`F(u)=0`, :math:`F=(F_0,F_1)`, :math:`u=(u_0,u_1)`, and compute the Jacobian :math:`J_{i,j}=\partial F_i/\partial u_j`. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:vib:geometric:mean: Problem 5.5: Solve nonlinear equations arising from a vibration ODE ------------------------------------------------------------------- Consider a nonlinear vibration problem .. _Eq:_auto252: .. math:: \tag{603} mu^{\prime\prime} + bu^{\prime}|u^{\prime}| + 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^{\prime\prime}` is mass times acceleration, and :math:`bu^{\prime}|u^{\prime}|` models water or air drag. **a)** Rewrite the equation for :math:`u` as a system of two first-order ODEs, and discretize this system by a Crank-Nicolson (centered difference) method. With :math:`v=u^\prime`, we get a nonlinear term :math:`v^{n+\frac{1}{2}}|v^{n+\frac{1}{2}}|`. Use a geometric average for :math:`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. .. 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 Filename: ``nonlin_vib``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:products:arith:mean: Exercise 5.6: Find the truncation error of arithmetic mean of products ---------------------------------------------------------------------- In the section :ref:`nonlin:pdelevel:Picard:CN` we introduce alternative arithmetic means of a product. Say the product is :math:`P(t)Q(t)` evaluated at :math:`t=t_{n+\frac{1}{2}}`. The exact value is .. math:: [PQ]^{n+\frac{1}{2}} = P^{n+\frac{1}{2}}Q^{n+\frac{1}{2}} There are two obvious candidates for evaluating :math:`[PQ]^{n+\frac{1}{2}}` as a mean of values of :math:`P` and :math:`Q` at :math:`t_n` and :math:`t_{n+1}`. Either we can take the arithmetic mean of each factor :math:`P` and :math:`Q`, .. _Eq:nonlin:exer:products:arith:mean:f: .. math:: \tag{604} [PQ]^{n+\frac{1}{2}} \approx \frac{1}{2}(P^n + P^{n+1})\frac{1}{2}(Q^n + Q^{n+1}), or we can take the arithmetic mean of the product :math:`PQ`: .. _Eq:nonlin:exer:products:arith:mean:p: .. math:: \tag{605} [PQ]^{n+\frac{1}{2}} \approx \frac{1}{2}(P^nQ^n + P^{n+1}Q^{n+1}){\thinspace .} The arithmetic average of :math:`P(t_{n+\frac{1}{2}})` is :math:`{\mathcal{O}(\Delta t^2)}`: .. math:: P(t_{n+\frac{1}{2}}) = \frac{1}{2}(P^n + P^{n+1}) +{\mathcal{O}(\Delta t^2)}{\thinspace .} A fundamental question is whether :ref:`(604) ` and :ref:`(605) ` have different orders of accuracy in :math:`\Delta t = t_{n+1}-t_n`. To investigate this question, expand quantities at :math:`t_{n+1}` and :math:`t_n` in Taylor series around :math:`t_{n+\frac{1}{2}}`, and subtract the true value :math:`[PQ]^{n+\frac{1}{2}}` from the approximations :ref:`(604) ` and :ref:`(605) ` to see what the order of the error terms are. .. --- begin hint in exercise --- **Hint.** You may explore ``sympy`` for carrying out the tedious calculations. A general Taylor series expansion of :math:`P(t+\frac{1}{2}\Delta t)` around :math:`t` involving just a general function :math:`P(t)` can be created as follows: .. code-block:: python >>> 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, :math:`\frac{1}{2}(P(-\frac{1}{2}\Delta t) + P(-\frac{1}{2}\Delta t))` for :math:`t=0` is then .. code-block:: python >>> 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 :ref:`(604) ` and :ref:`(605) ` for :math:`n=0`. (Choosing :math:`n=0` is necessary for not making the expressions too complicated for ``sympy``, but there is of course no lack of generality by using :math:`n=0` rather than an arbitrary :math:`n` - the main point is the product and addition of Taylor series.) .. --- end hint in exercise --- Filename: ``product_arith_mean``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:Newton:linear: Problem 5.7: 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: ``Newton_linear``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:1pu2:fem: Problem 5.8: Discretize a 1D problem with a nonlinear coefficient ----------------------------------------------------------------- We consider the problem .. _Eq:nonlin:exer:1D:1pu2:fem:pde: .. math:: \tag{606} ((1 + u^2)u^{\prime})^{\prime} = 1,\quad x\in (0,1),\quad u(0)=u(1)=0{\thinspace .} Discretize :ref:`(606) ` by a centered finite difference method on a uniform mesh. Filename: ``nonlin_1D_coeff_discretize``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:1pu2:PicardNewton: Problem 5.9: Linearize a 1D problem with a nonlinear coefficient ---------------------------------------------------------------- We have a two-point boundary value problem .. _Eq:nonlin:exer:1D:1pu2:PicardNewton:pde: .. math:: \tag{607} ((1 + u^2)u^{\prime})^{\prime} = 1,\quad x\in (0,1),\quad u(0)=u(1)=0{\thinspace .} **a)** Construct a Picard iteration method for :ref:`(607) ` without discretizing in space. **b)** Apply Newton's method to :ref:`(607) ` without discretizing in space. **c)** Discretize :ref:`(607) ` by a centered finite difference scheme. Construct a Picard method for the resulting system of nonlinear algebraic equations. **d)** Discretize :ref:`(607) ` 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``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:fu:discretize:fd: Problem 5.10: Finite differences for the 1D Bratu problem --------------------------------------------------------- We address the so-called Bratu problem .. _Eq:nonlin:exer:1D:fu:discretize:fd:pde: .. math:: \tag{608} u^{\prime\prime} + \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 :ref:`(608) ` has an exact solution .. math:: {u_{\small\mbox{e}}}(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 :ref:`(608) ` 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 :ref:`(608) ` 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. **c)** Implement a solver that can compute :math:`u(x)` using Newton's method. Plot the error as a function of :math:`x` in each iteration. **d)** Investigate whether Newton's method gives second-order convergence by computing :math:`|| {u_{\small\mbox{e}}} - u||/||{u_{\small\mbox{e}}} - u^{-}||^2` in each iteration, where :math:`u` is solution in the current iteration and :math:`u^{-}` is the solution in the previous iteration. Filename: ``nonlin_1D_Bratu_fd``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:1D:heat:nonlinear:fdm: Problem 5.11: Discretize a nonlinear 1D heat conduction PDE by finite differences --------------------------------------------------------------------------------- We address the 1D heat conduction PDE .. math:: \varrho c(T) T_t = (k(T)T_x)_x, for :math:`x\in [0,L]`, 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. :math:`T(x,0)=I(x)`, and ends are subject to a cooling law: .. math:: k(T)T_x|_{x=0} = h(T)(T-T_s),\quad -k(T)T_x|_{x=L}=h(T)(T-T_s), where :math:`h(T)` is a heat transfer coefficient and :math:`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``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:grad:pow:term: Problem 5.12: 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. The expression :math:`|\nabla u|` is defined as the Euclidean norm of a vector: :math:`|\nabla u|^2 = \nabla u \cdot \nabla u`. 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 .} .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``nonlin_differentiate``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:2D:heat:nonlinear:fd: Exercise 5.13: Crank-Nicolson for a nonlinear 3D diffusion equation ------------------------------------------------------------------- Redo the section :ref:`nonlin:alglevel:dD:fd` when a Crank-Nicolson scheme is used to discretize the equations in time and the problem is formulated for three spatial dimensions. .. --- begin hint in exercise --- **Hint.** Express the Jacobian as :math:`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 :math:`J_{i,j,k,r,s,t}` is very sparse: :math:`J_{i,j,k,r,s,t}\neq 0` only for :math:`r=i\pm i`, :math:`s=j\pm 1`, and :math:`t=k\pm 1` as well as :math:`r=i`, :math:`s=j`, and :math:`t=k`. .. --- end hint in exercise --- Filename: ``nonlin_heat_FD_CN_2D``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:sparsity:Jacobian: Problem 5.14: 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 equation at a point :math:`(i,j)` in 2D or :math:`(i,j,k)` in 3D, and identify the nonzero entries of the Jacobian that can arise from such a type of difference equation. .. --- end hint in exercise --- Filename: ``nonlin_sparsity_Jacobian``. .. --- end exercise --- .. --- begin exercise --- .. _nonlin:exer:continuation:1DnNflow: Problem 5.15: 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^{\prime}(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 ---