# 3.7 Nonlinear problems¶

We want to solve a nonlinear PDE.

## A simple scalar PDE¶

We consider the simple PDE

$- \Delta u + 3 u^3 = 1 \text{ in } \Omega$

on the unit square $$\Omega = (0,1)^2$$.

We note that this PDE can also be formulated as a nonlinear minimization problem (cf. 3.8).

:

from ngsolve import *
from netgen.geom2d import unit_square
from ngsolve.webgui import Draw
mesh = Mesh (unit_square.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

$A(u,v) = \int_{\Omega} \nabla u \nabla v + 3 u^3 v - 1 v ~ dx \quad ( = 0 ~ \forall~v \in H^1_0)$
:

V = H1(mesh, order=3, dirichlet=[1,2,3,4])
u,v = V.TnT()
a = BilinearForm(V)


### Newton’s method¶

We use Newton’s method and make the loop:

• 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$$.

:

def SimpleNewtonSolve(gfu,a,tol=1e-13,maxits=25):
res = gfu.vec.CreateVector()
du = gfu.vec.CreateVector()
fes = gfu.space
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

#stopping criteria
stopcritval = sqrt(abs(InnerProduct(du,res)))
print ("<A u",it,", A u",it,">_{-1}^0.5 = ", stopcritval)
if stopcritval < tol:
break

:

gfu = GridFunction(V)
Draw(gfu,mesh,"u")
SimpleNewtonSolve(gfu,a)

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  0.18743829125307696
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  9.417800751716788e-05
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  8.541508292330884e-11
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  4.019346184174276e-17


There are also some solvers shipped with NGSolve now:

:

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='umfpack', 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.


:

gfu.vec[:]=0
Newton(a,gfu,freedofs=gfu.space.FreeDofs(),maxit=100,maxerr=1e-11,inverse="umfpack",dampfactor=1,printing=True)

Newton iteration  0
err =  0.18743829125307693
Newton iteration  1
err =  9.417800751713952e-05
Newton iteration  2
err =  8.541507396334433e-11
Newton iteration  3
err =  3.5543860698953376e-17

:

(0, 4)


## A trivial problem:¶

$5 u^2 = 1, \qquad u \in \mathbb{R}.$
:

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)

print("\nscalar solution", gfu.vec, "(exact: ", sqrt(0.2), ")")

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  1.2649110640673513
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.32659863237109055
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  0.04114755998989123
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  0.0008574269268691781
Iteration   4  <A u 4 , A u 4 >_{-1}^0.5 =  3.883274522609999e-07
Iteration   5  <A u 5 , A u 5 >_{-1}^0.5 =  7.979879233426057e-14

scalar solution 0.4472135954999579 (exact:  0.4472135954999579 )


## Another example: Stationary Navier-Stokes:¶

Find $$\mathbf{u} \in \mathbf{V}$$, $$p \in Q$$, $$\lambda \in \mathbb{R}$$ so that \begin{align} \int_{\Omega} \nu \nabla \mathbf{u} : \nabla \mathbf{v} + (\mathbf{u} \cdot \nabla) \mathbf{u} \cdot \mathbf{v}& - \int_{\Omega} \operatorname{div}(\mathbf{v}) p & &= \int \mathbf{f} \cdot \mathbf{v} && \forall \mathbf{v} \in \mathbf{V}, \\ - \int_{\Omega} \operatorname{div}(\mathbf{u}) q & & + \int_{\Omega} \lambda q &= 0 && \forall q \in Q, \\ & \int_{\Omega} \mu p & &= 0 && \forall \mu \in \mathbb{R}. \end{align}

:

mesh = Mesh (unit_square.GenerateMesh(maxh=0.05)); nu = Parameter(1)
V = VectorH1(mesh,order=3,dirichlet="bottom|right|top|left")
Q = H1(mesh,order=2);
N = NumberSpace(mesh);
X = V*Q*N
(u,p,lam), (v,q,mu) = X.TnT()
a = BilinearForm(X)
-div(u)*q-div(v)*p-lam*q-mu*p)*dx

:

gfu = GridFunction(X)
gfu.components.Set(CoefficientFunction((4*x*(1-x),0)),
definedon=mesh.Boundaries("top"))

:

SimpleNewtonSolve(gfu,a)
scenep = Draw(gfu.components,mesh,"p")
sceneu = Draw(gfu.components,mesh,"u")

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  2.7970749709609817
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.007939103696589463
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  6.512315396348492e-08
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  8.92989109016202e-16

:

nu.Set(0.01)
SimpleNewtonSolve(gfu,a)
sceneu.Redraw()
scenep.Redraw()

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  0.08811928591680591
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.008271027344739551
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  0.00010638956077082668
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  1.9401670843845306e-08
Iteration   4  <A u 4 , A u 4 >_{-1}^0.5 =  5.786473673734125e-16

:

nu.Set(0.001)
SimpleNewtonSolve(gfu,a)
sceneu.Redraw()
scenep.Redraw()

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  0.07337564494746804
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.05991419000703266
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  0.036275704108326874
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  0.031863369330183786
Iteration   4  <A u 4 , A u 4 >_{-1}^0.5 =  0.030208729410974255
Iteration   5  <A u 5 , A u 5 >_{-1}^0.5 =  0.024437032676887488
Iteration   6  <A u 6 , A u 6 >_{-1}^0.5 =  0.05105777725377034
Iteration   7  <A u 7 , A u 7 >_{-1}^0.5 =  0.10764886843571098
Iteration   8  <A u 8 , A u 8 >_{-1}^0.5 =  0.2550606781247994
Iteration   9  <A u 9 , A u 9 >_{-1}^0.5 =  0.4025326927154164
Iteration  10  <A u 10 , A u 10 >_{-1}^0.5 =  0.544680529442505
Iteration  11  <A u 11 , A u 11 >_{-1}^0.5 =  0.9286504038561297
Iteration  12  <A u 12 , A u 12 >_{-1}^0.5 =  7.115744662274641
Iteration  13  <A u 13 , A u 13 >_{-1}^0.5 =  50.20546195200112
Iteration  14  <A u 14 , A u 14 >_{-1}^0.5 =  31.943951750267242
Iteration  15  <A u 15 , A u 15 >_{-1}^0.5 =  256.4716270013097
Iteration  16  <A u 16 , A u 16 >_{-1}^0.5 =  78.06232795008292
Iteration  17  <A u 17 , A u 17 >_{-1}^0.5 =  25.783039614789075
Iteration  18  <A u 18 , A u 18 >_{-1}^0.5 =  409.3711287752754
Iteration  19  <A u 19 , A u 19 >_{-1}^0.5 =  2636.9491272880073
Iteration  20  <A u 20 , A u 20 >_{-1}^0.5 =  1504.1151501285383
Iteration  21  <A u 21 , A u 21 >_{-1}^0.5 =  1188.8074227234345
Iteration  22  <A u 22 , A u 22 >_{-1}^0.5 =  3215.624363281052
Iteration  23  <A u 23 , A u 23 >_{-1}^0.5 =  7663.189625717757
Iteration  24  <A u 24 , A u 24 >_{-1}^0.5 =  2576.0342401692915

:

nu.Set(0.001)
gfu.components.Set(CoefficientFunction((4*x*(1-x),0)),definedon=mesh.Boundaries("top"))
Newton(a,gfu,maxit=20,dampfactor=0.1)
sceneu.Redraw()
scenep.Redraw()

Newton iteration  0
err =  2.7739042910681238
Newton iteration  1
err =  2.496506225593245
Newton iteration  2
err =  1.997288467071201
Newton iteration  3
err =  1.3983527685940964
Newton iteration  4
err =  0.8396015387654532
Newton iteration  5
err =  0.4216112662057235
Newton iteration  6
err =  0.17108085772525564
Newton iteration  7
err =  0.051789307422249006
Newton iteration  8
err =  0.010849749736353269
Newton iteration  9
err =  0.0012114726993057965
Newton iteration  10
err =  1.258495602389179e-05
Newton iteration  11
err =  6.675805454868902e-09
Newton iteration  12
err =  1.5702971755330574e-15

[ ]: