Nonlinear 1d PDE

Linear Algebra
1

Dear all,
I want to solve the 1d nonlinear differential equation \frac {d^2 u} {dz^2} =-\theta(z)u+\frac {1} {u^3} +\frac {Q} {2u}, where \theta(z) is the step function and Q is set to 10. The code written with dolfinx is shown below:

import numpy as np

from mpi4py import MPI

from petsc4py.PETSc import ScalarType, Options

from dolfinx import mesh, fem, io, nls

from ufl import TestFunction, inner, grad, dx

domain = mesh.create_interval(MPI.COMM_WORLD, 20, [0, 1])

V = fem.FunctionSpace(domain, ("CG", 1))

u = fem.Function(V)

v = TestFunction(V)

theta = fem.Function(V)

cells_0 = mesh.locate_entities(domain, domain.topology.dim, lambda x: x[0]<=0.5)

cells_1 = mesh.locate_entities(domain, domain.topology.dim, lambda x: x[0]>=0.5)

theta.x.array[cells_0] = np.full_like(cells_0, 1, dtype=ScalarType)

theta.x.array[cells_1] = np.full_like(cells_1, 0, dtype=ScalarType)

F = inner(u**3*grad(u), grad(v))*dx - inner(theta*u**4, v)*dx + inner(fem.Constant(domain, ScalarType(1)), v)*dx + inner(10*u**2, v)*dx

problem = fem.petsc.NonlinearProblem(F, u)

solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)

solver.convergence_criterion = "incremental"

solver.rtol = 1e-4

solver.report = True

n, converged = solver.solve(u)

assert(converged)

print(f"Number of interations: {n:d}")

However, an error is reported:

RuntimeError: Newton solver did not converge because maximum number of iterations reached.

I’m really appreciated if anyone can help me.