Page 40 - Process Modelling and Simulation With Finite Element Methods

FEMLAB and the Basics of Numerical Analysis 21

Note in proof I got into the habit of using the coefficient mode and numeric
Jacobian as it mirrors my style of FEM coding - the coefficient mode does not
include all the potential nonlinearity in the coefficients in assembling the
stiffness matrix, but by computing the Jacobian numerically, it is all included in
the iterative scheme. For small problems, you will see no performance
degradation. On larger problems, it is wiser to use the PDE general mode and
the exact Jacobian option, which assembles the full nonlinear contributions to the
Jacobian analytically. If I had to write this chapter again, all the coefficient
modes would disappear.
During the Solve step, the report window shows a runout of several columns
per iteration; particularly important is the error estimate ErrEst, which for
iteration 4 is about 10.’ which is smaller in magnitude than the default
tolerance set at in the Nonlinear tab of the Solver Parameters dialog
box. Click anywhere on the grid and the report window will “Value of u(u) at
0.456: -2.73205.” The analytically determined root nearest to this is -I-&,
showing the numerical solution in good agreement. According to the structure of
the quadratic formula of algebra, clearly another root is -I+&, and by
inspection, the third root is 1. Returning to the subdomain settings, set the initial
guess to u(tO)=-0.5 and FEMLAB converges to u=0.732051, again a good
approximation. u(tO)=l.2 as an initial guess converges to u=l.
This exercise clarifies two features of nonlinear solvers and problems - (i)
nonlinear problems can have multiple solutions; (ii) the initial guess is key to
convergence to a particular solution. With Newton’s method, it is usually the
case that convergence is to the nearest solution, but overshoots in highly
nonlinear problems may override this stereotype. These features persist in
higher dimensional solution spaces and with spatial-temporal dependence.
The MATLAB model m-file ro0tfinder.m contains all the MATLAB source
code with FEMLAB extensions to reproduce the current state of the FEMLAB
GUI. This file is available from the website http://eyrie.shef.ac.uk/femlab. Just
pull down the file menu, select Open model m-file, and use the Open file dialog
window to locate it. You can rapidly place your nonlinear function in the
Subdomain settings, specify an initial guess, and use the stationary nonlinear
solver to converge to a solution. But what if your function does not have a linear
component to put on the LHS of (1.3)? For instance, tanh(u) - u2 + 5 = 0
results in a singular stiffness matrix when FEMLAB assembles the LHS of (1.3).
The suggestion is to set the coefficient of the second derivative of u, c=l in the
Subdomain settings. Coupled with the Neumann boundary conditions, this
artificial diffusion cannot change the fact that the solution must be constant over
the single element, yet it prevents the stiffness matrix from becoming singular.

Root Finding in General Mode
The difficulty with a singular stiffness matrix assembly for tanh(u) - u2 + 5 = 0 can
be averted by using General Mode, which solves