This page was generated from 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).
[1]:
from ngsolve import *
import netgen.gui
from netgen.geom2d import unit_square
mesh = Mesh (unit_square.GenerateMesh(maxh=0.3))
Generate Mesh from spline geometry
Boundary mesh done, np = 12
CalcLocalH: 12 Points 0 Elements 0 Surface Elements
Meshing domain 1 / 1
load internal triangle rules
Surface meshing done
Edgeswapping, topological
Smoothing
Split improve
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Split improve
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Split improve
Combine improve
Smoothing
Update mesh topology
Update clusters
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)\]
[2]:
V = H1(mesh, order=3, dirichlet=[1,2,3,4])
u,v = V.TnT()
a = BilinearForm(V)
a += (grad(u) * grad(v) + 3*u**3*v- 1 * v)*dx
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\).
[3]:
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
[4]:
gfu = GridFunction(V)
Draw(gfu,mesh,"u")
SimpleNewtonSolve(gfu,a)
Iteration 0 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 0 , A u 0 >_{-1}^0.5 = 0.187435321076685
Iteration 1 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 1 , A u 1 >_{-1}^0.5 = 9.416472885501191e-05
Iteration 2 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 2 , A u 2 >_{-1}^0.5 = 8.538473326331901e-11
Iteration 3 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 3 , A u 3 >_{-1}^0.5 = 3.831924934146655e-17
There are also some solvers shipped with NGSolve now:
[5]:
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.
[6]:
gfu.vec[:]=0
Newton(a,gfu,freedofs=gfu.space.FreeDofs(),maxit=100,maxerr=1e-11,inverse="umfpack",dampfactor=1,printing=True)
Newton iteration 0
Assemble linearization
assemble VOL element 24/24
call umfpack ... done
err = 0.18743532107668498
Newton iteration 1
Assemble linearization
assemble VOL element 24/24
call umfpack ... done
err = 9.416472885497372e-05
Newton iteration 2
Assemble linearization
assemble VOL element 24/24
call umfpack ... done
err = 8.538474393723481e-11
Newton iteration 3
Assemble linearization
assemble VOL element 24/24
call umfpack ... done
err = 4.291958722351373e-17
[6]:
(0, 4)
A trivial problem:¶
\[5 u^2 = 1, \qquad u \in \mathbb{R}.\]
[7]:
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[0], "(exact: ", sqrt(0.2), ")")
Iteration 0 Assemble linearization
assemble VOL element 24/24
call pardiso ...<A u done
0 , A u 0 >_{-1}^0.5 = 1.2649110640673515
Iteration 1 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 1 , A u 1 >_{-1}^0.5 = 0.32659863237109055
Iteration 2 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 2 , A u 2 >_{-1}^0.5 = 0.04114755998989123
Iteration 3 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 3 , A u 3 >_{-1}^0.5 = 0.0008574269268691781
Iteration 4 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 4 , A u 4 >_{-1}^0.5 = 3.8832745226099997e-07
Iteration 5 Assemble linearization
assemble VOL element 24/24
call pardiso ... done
<A u 5 , A u 5 >_{-1}^0.5 = 7.979879233426059e-14
scalar solution 0.4472135954999579 (exact: 0.4472135954999579 )