This page was generated from unit-3.3-nonlinear/nonlinear.ipynb.
3.3 Nonlinear problems¶
In this unit we turn our attention to nonlinear PDE problems and the tools that NGSolve provides to simplify it.
A simple scalar PDE¶
Let us start with a simple PDE with a nonlinearity:
on the unit square \(\Omega = (0,1)^2\).
We note that this PDE can also be formulated as a nonlinear minimization problem (cf. 3.4).
[1]:
# define geometry and generate mesh
from ngsolve import *
from ngsolve.webgui import *
from netgen.occ import *
shape = Rectangle(1,1).Face()
shape.edges.Min(X).name="left"
shape.edges.Max(X).name="right"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(Y).name="top"
geom = OCCGeometry(shape, dim=2)
mesh = Mesh(geom.GenerateMesh(maxh=0.3))
In NGSolve we can solve the PDE conveniently using the linearization feature of SymbolicBFI.
The BilinearForm (which is not bilinear!) needed in the weak formulation is
[2]:
V = H1(mesh, order=3, dirichlet=[1,2,3,4])
u,v = V.TnT()
a = BilinearForm(V)
a += (grad(u) * grad(v) + 1/3*u**3*v- 10 * v)*dx
Newton’s method¶
In the preparation of the formulation of Newton’s method we identify the semilinear form \(A(\cdot,\cdot)\) with the corresponding operator \(A: V \to V'\) and write (assuming sufficient regularity)
\(A(u)\) for the linear form \([A(u)](v) = A(u,v)\) and
the derivative of \(A\) at an evaluation point \(u\) is the bilinear form \(\delta A(u)\) with
\[[\delta A(u)] (w,v) = \lim_{h \to 0} \frac{ A(u+h\cdot w,v) - A(u,v)}{h}\]
You could compute these linear and bilinear forms manually. However, NGSolve provides functions for both operations:
The linear form \(A(u)\) for a given (vector) u is obtained from
a.Apply(vec,res)(resulting in a vectorresrepresenting the linear form)The bilinear form (represented by the corresponding matrix) is stored in
a.matafter a call ofa.AssembleLinearization(vec)
Under the hood it uses its functionality to derive CoefficientFunctions symbolically.
For example, NGSolve derives the bilinear form integrand $ \frac13 u^3 v $ w.r.t. \(u\) at \(\tilde{u}\) in direction \(w\) resulting in the integrand \(\tilde{u}^2 w v\).
This allows to form the corresponding bilinear form integrals automatically for you.
To obtain a Newton algorithm we hence only need to translate the following pseudo-code formulation of Newton’s method to NGSolve. The pseudo code is:
NewtonSolve(Pseudo code):
Given an initial guess \(u^0\)
loop over \(i=0,..\) until convergence:
Compute linearization: \(A (u^i) + \delta A(u^i) \Delta u^{i} = 0\):
\(f^i = A (u^i)\)
\(B^i = \delta A(u^i)\)
Solve \(B^i \Delta u^i = -f^i\)
Update \(u^{i+1} = u^i + \Delta u^{i}\)
Evaluate stopping criteria
As a stopping criteria we take \(\langle A u^i,\Delta u^i \rangle = \langle A u^i, A u^i \rangle_{(B^i)^{-1}}< \varepsilon\).
Now, here comes the same thing in NGSolve syntax:
[3]:
def SimpleNewtonSolve(gfu,a,tol=1e-13,maxits=10, callback=lambda gfu: None):
res = gfu.vec.CreateVector()
du = gfu.vec.CreateVector()
fes = gfu.space
callback(gfu)
for it in range(maxits):
print ("Iteration {:3} ".format(it),end="")
a.Apply(gfu.vec, res)
a.AssembleLinearization(gfu.vec)
du.data = a.mat.Inverse(fes.FreeDofs()) * res
gfu.vec.data -= du
callback(gfu)
#stopping criteria
stopcritval = sqrt(abs(InnerProduct(du,res)))
print ("<A u",it,", A u",it,">_{-1}^0.5 = ", stopcritval)
if stopcritval < tol:
break
Let’s apply this to the previous PDE problem:
[4]:
gfu = GridFunction(V)
gfu.Set((x*(1-x))**4*(y*(1-y))**4) # initial guess
gfu_it = GridFunction(gfu.space,multidim=0)
cb = lambda gfu : gfu_it.AddMultiDimComponent(gfu.vec) # store current state
SimpleNewtonSolve(gfu, a, callback = cb)
Iteration 0 <A u 0 , A u 0 >_{-1}^0.5 = 1.8743634219605203
Iteration 1 <A u 1 , A u 1 >_{-1}^0.5 = 0.010378681410246811
Iteration 2 <A u 2 , A u 2 >_{-1}^0.5 = 1.122416328991341e-06
Iteration 3 <A u 3 , A u 3 >_{-1}^0.5 = 1.3180540078433946e-14
[5]:
Draw(gfu,mesh,"u", deformation = True)
[5]:
BaseWebGuiScene
[6]:
Draw(gfu_it,mesh,"u", deformation = True)
[6]:
BaseWebGuiScene
Use the multidim-Slider to inspect the results after the iterations.
There are also some solvers shipped with NGSolve now. The ngsolve.solvers.Newton method allows you to use static condensation and is overall more refined (e.g. with damping options). For this tutorial we will mostly stay with the simple hand-crafted SimpleNewtonSolve. Here is the alternative demonstrated once, nevertheless:
[7]:
from ngsolve.solvers import *
help(Newton)
Help on function Newton in module ngsolve.nonlinearsolvers:
Newton(a, u, freedofs=None, maxit=100, maxerr=1e-11, inverse='', dirichletvalues=None, dampfactor=1, printing=True, callback=None)
Newton's method for solving non-linear problems of the form A(u)=0.
Parameters
----------
a : BilinearForm
The BilinearForm of the non-linear variational problem. It does not have to be assembled.
u : GridFunction
The GridFunction where the solution is saved. The values are used as initial guess for Newton's method.
freedofs : BitArray
The FreeDofs on which the assembled matrix is inverted. If argument is 'None' then the FreeDofs of the underlying FESpace is used.
maxit : int
Number of maximal iteration for Newton. If the maximal number is reached before the maximal error Newton might no converge and a warning is displayed.
maxerr : float
The maximal error which Newton should reach before it stops. The error is computed by the square root of the inner product of the residuum and the correction.
inverse : string
A string of the sparse direct solver which should be solved for inverting the assembled Newton matrix.
dampfactor : float
Set the damping factor for Newton's method. If dampfactor is 1 then no damping is done. If value is < 1 then the damping is done by the formula 'min(1,dampfactor*numit)' for the correction, where 'numit' denotes the Newton iteration.
printing : bool
Set if Newton's method should print informations about the actual iteration like the error.
Returns
-------
(int, int)
List of two integers. The first one is 0 if Newton's method did converge, -1 otherwise. The second one gives the number of Newton iterations needed.
We call this Newton method as well here.
[8]:
gfu.Set((x*(1-x))**4*(y*(1-y))**4) # initial guess
Newton(a,gfu,freedofs=gfu.space.FreeDofs(),maxit=100,maxerr=1e-11,inverse="umfpack",dampfactor=1,printing=True)
Newton iteration 0
err = 1.87436342196052
Newton iteration 1
err = 0.010378681410246426
Newton iteration 2
err = 1.1224163290972609e-06
Newton iteration 3
err = 1.3421225797380163e-14
[8]:
(0, 4)
A trivial problem:¶
As a second simple problem, let us consider a trivial scalar problem:
We chose this problem and put this in the (somewhat artificially) previous PDE framework and solve with Newton’s method in a setting that you could easily follow with pen and paper:
[9]:
V = NumberSpace(mesh)
u,v = V.TnT()
a = BilinearForm(V)
a += ( 5*u*u*v - 1 * v)*dx
gfu = GridFunction(V)
gfu.vec[:] = 1
SimpleNewtonSolve(gfu,a, callback = lambda gfu : print(f"u^k = {gfu.vec[0]}, u^k**2 = {gfu.vec[0]**2}"))
print(f"\nscalar solution, {gfu.vec[0]}, exact: {sqrt(0.2)}, error: {abs(sqrt(0.2)-gfu.vec[0])}")
u^k = 1.0, u^k**2 = 1.0
Iteration 0 u^k = 0.6, u^k**2 = 0.36
<A u 0 , A u 0 >_{-1}^0.5 = 1.264911064067352
Iteration 1 u^k = 0.4666666666666667, u^k**2 = 0.2177777777777778
<A u 1 , A u 1 >_{-1}^0.5 = 0.3265986323710903
Iteration 2 u^k = 0.4476190476190476, u^k**2 = 0.20036281179138318
<A u 2 , A u 2 >_{-1}^0.5 = 0.04114755998989124
Iteration 3 u^k = 0.4472137791286727, u^k**2 = 0.20000016424254927
<A u 3 , A u 3 >_{-1}^0.5 = 0.000857426926869178
Iteration 4 u^k = 0.4472135954999957, u^k**2 = 0.20000000000003376
<A u 4 , A u 4 >_{-1}^0.5 = 3.8832745226099997e-07
Iteration 5 u^k = 0.4472135954999579, u^k**2 = 0.19999999999999998
<A u 5 , A u 5 >_{-1}^0.5 = 7.979879233426059e-14
scalar solution, 0.4472135954999579, exact: 0.4472135954999579, error: 0.0