A Semi-Linear Problem
In the following we discuss a non-linear example. The so-called standard, semilinear PDE given by
\[\begin{aligned} -\Delta y + y^3 = f &\quad \text{in}\; \Omega\\ y = 0 &\quad \text{on}\; \Gamma. \end{aligned}\]
The variational form reads as follows: Find $y \in H_0^1(\Omega)$ such that
\[\int_\Omega \nabla y \cdot \nabla v\, dx + \int_\Omega y^3 v\, dx = \int_\Omega f v\, dx\]
for all $v \in H_0^1(\Omega)$.
Note that this is a non-linear equation in $y$. In order to apply Newton's method we have to differentiate with respect to $y$.
The linearization evaluated in $y$ in direction $w$ is given by
\[\int_\Omega \nabla w \cdot \nabla v\, dx + 3 \int_\Omega y^2 w v\, dx.\]
Consequently, one step of Newton's method is given by the solution of the variational problem: Find $w \in H^1_0(\Omega)$ such that
\[\int_\Omega \nabla w \cdot \nabla v\, dx + 3 \int_\Omega y^2 w v\, dx = -\int_\Omega \nabla y \cdot \nabla v\, dx - \int_\Omega y^3 v\, dx + \int_\Omega f v\, dx\]
for all $v \in H_0^1(\Omega)$.
We thus introduce the function assemble_cubicterm and assemble_cubicderivativematrix in order to assemble the source term involving $y^3$ and the matrix with $y^2$, respectively.
Most of the code is similar to the Poisson example. We thus illuminate what is required additionally. The function semilinear
evaluates Newton's method for the problem described above with a tolerance tol
:
using MinFEM, LinearAlgebra
function semilinear(mesh::Mesh, L::AbstractMatrix, M::AbstractMatrix,
s::AbstractVector, boundaryIndices::Set{Int64};
tol::Float64=1e-10, maxIterations::Int64=10)
y = zeros(mesh.nnodes)
pde = PDESystem(A=L, b=M*s, bc=zeros(mesh.nnodes), DI=boundaryIndices)
for i = 1:maxIterations
pde.A = L + assemble_cubicderivativematrix(mesh, y)
pde.b = -L*y + M*s - assemble_cubicterm(mesh, y)
refresh!(pde)
solve!(pde)
y += pde.state
res = norm(pde.state)
if res < tol
println("It. $i: $res < $tol")
println("Semi-linear routine converged.")
break
else
println("It. $i: $res ≥ $tol")
if i == maxIterations
println("Semi-linear routine failed.")
println("Maximum number of iterations reached.")
end
end
end
return y
end
The rest of the code is then similar to the Poisson problem:
mesh = import_mesh("../meshes/Zshaped.msh")
L = assemble_laplacian(mesh)
M = assemble_massmatrix(mesh)
f(x) = 3*2*pi^2*sin(x[1]*pi)*sin(x[2]*pi) + (3*sin(x[1]*pi)*sin(x[2]*pi))^3
s = evaluate_mesh_function(mesh, f)
boundary = select_boundaries(mesh)
boundaryNodes = extract_nodes(boundary)
y = semilinear(mesh, L, M, s, boundaryNodes)
write_to_vtk([y, s], mesh, ["y","s"], "semilinear")
In Paraview, the visualization should then look similar to the following: