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 *
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)\]
[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 <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:
[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
err = 0.18743829125307693
Newton iteration 1
err = 9.417800751713952e-05
Newton iteration 2
err = 8.541507396334433e-11
Newton iteration 3
err = 3.5543860698953376e-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 <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 )