We want to solve a nonlinear PDE.
import netgen.gui
%gui tk
import tkinter
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 + 5 u^2 v - 1 v ~ dx \quad ( = 0 ~ \forall~v \in H^1_0) $$
V = H1(mesh, order=3, dirichlet=[1,2,3,4])
u = V.TrialFunction()
v = V.TestFunction()
a = BilinearForm(V)
a += SymbolicBFI( grad(u) * grad(v) + 5*u*u*v- 1 * v)
We use Newton's method and make the loop:
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$.
u = GridFunction(V)
Draw(u,mesh,"u")
res = u.vec.CreateVector()
du = u.vec.CreateVector()
for it in range(25):
print ("Iteration",it)
a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)
#a.Assemble()
du.data = a.mat.Inverse(V.FreeDofs()) * res
u.vec.data -= du
#stopping criteria
stopcritval = sqrt(abs(InnerProduct(du,res)))
print ("<A u",it,", A u",it,">_{-1}^0.5 = ",
stopcritval)
if stopcritval < 1e-13:
break