This page was generated from jupyter-files/unit-3.7-nonlinear/nonlinear.ipynb.

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).

import netgen.gui
%gui tk
from netgen.geom2d import unit_square
from ngsolve import *

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)
a += SymbolicBFI( grad(u) * grad(v) + 5*u*u*v- 1 * 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.18745474302207424
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.002571436669781315
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  5.326881678156768e-07
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  2.3069719054272302e-14

A trivial problem:

\[5 u^2 = 1, \qquad u \in \mathbb{R}.\]

V = NumberSpace(mesh)
u,v = V.TnT()
a = BilinearForm(V)
a += SymbolicBFI( 5*u*u*v - 1 * v)
gfu = GridFunction(V)
gfu.vec[:] = 1
SimpleNewtonSolve(gfu,a)

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

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  1.2649110640673518
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.3265986323710904
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  0.041147559989891225
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.8832745226099987e-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))
V = VectorH1(mesh,order=3, dirichlet="bottom|right|top|left")
Q = H1(mesh,order=2)
N = NumberSpace(mesh)
X = FESpace([V,Q,N])
(u,p,lam), (v,q,mu) = X.TnT()
a = BilinearForm(X)
nu = Parameter(1)
# Note: grad(u) is a matrix where each row is the gradient of one of the components
a += SymbolicBFI(nu*InnerProduct(grad(u),grad(v))+InnerProduct(grad(u)*u,v)-div(u)*q-div(v)*p-lam*q-mu*p)
gfu = GridFunction(X)
gfu.components[0].Set(CoefficientFunction((4*x*(1-x),0)),definedon=mesh.Boundaries("top"))

SimpleNewtonSolve(gfu,a)
Draw(gfu.components[1],mesh,"p")
Draw(gfu.components[0],mesh,"u")

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  2.8037383397141853
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.00795105147804976
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  6.527752846290048e-08
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  9.248737238644841e-16

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

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  0.08811927338814224
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.008271022172200517
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  0.00010638943511847134
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  1.9401631614310508e-08
Iteration   4  <A u 4 , A u 4 >_{-1}^0.5 =  5.738505921143496e-16