# Newton solver does not converge due to non-linear term

**URL:** <https://fenicsproject.discourse.group/t/newton-solver-does-not-converge-due-to-non-linear-term/4176>\
**Category:** Uncategorized\
**Created:** [September 26, 2020, 7:32pm UTC](https://fenicsproject.discourse.group/t/newton-solver-does-not-converge-due-to-non-linear-term/4176 "2020-09-26T19:32:54Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![spoulain](https://avatars.discourse-cdn.com/v4/letter/s/0ea827/32.png) [@spoulain](https://fenicsproject.discourse.group/u/spoulain)\
**Post date:** [September 26, 2020, 7:32pm UTC](https://fenicsproject.discourse.group/t/newton-solver-does-not-converge-due-to-non-linear-term/4176/1 "2020-09-26T19:32:54Z")

</div>

Hello,

I am very new to FEniCS and trying to solve the following time-dependent, non-linear, 1D equation:  
\frac{\partial h}{\partial t} + \frac{\partial}{\partial x}\left(h^3\frac{\partial^3h}{\partial x^3} + \frac 1h \frac{\partial h}{\partial x} \right)=0  
with boundary conditions:  
\frac{\partial h}{\partial x}=\frac{\partial^3 h}{\partial x^3}=0 ~~\text{at}~~ x=0,\frac{\lambda}{2}  
and initial condition:  
h(x,t=0)=\frac{10}{9} \left(1-\frac{1}{10}\cos(\frac{2\pi x}{\lambda})\right)

Introducing p=\frac{\partial^2 h}{\partial x^2} it is straightforward to obtain the variational formulation, with homogeneous Neumann boundary conditions for h and p. For now I only use a simple forward Euler scheme for the time.

I can’t do anything with that equation in FEniCS because of the term \frac{1}{h} \frac{\partial h}{\partial x}: at the first time-iteration, the newton algorithm does not converge and in fact the residual is already NaN at iteration 0. Everything seems fine when I remove it, or even replace the 1/h with some other functions of h.

I do not understand what could cause this. I do expect 1/h to be a issue eventually and lead to a singularity: h(x=0,t) \rightarrow 0 and I will have to deal with that later, but it should be OK at early times since the initial condition is nowhere close to 0.

I searched a bit in this forum but did not find tricks that could work, except using a random initial guess: I do not understand why but this seemingly improves things a bit: the Newton solver doesn’t converge, but at least the residual starts with a reasonable, non-NaN value (I tried to adjust the relaxation parameter, tolerances, … but couldn’t get it to converge).

Thanks for any help. My code is below

Note: This is a problem taken from the following reference  
Zhang, W.W. and Lister, J.R., 1999. Similarity solutions for van der Waals rupture of a thin film on a solid substrate. _Physics of Fluids_ , _11_ (9), pp.2454-2462.

```
from fenics import *
import numpy as np
from numpy.random import rand

t_max = 20
DT = 0.001
dt = Constant(DT)

lambd=10
h0=Expression('10./9. * (1-0.1*cos(2*3.1415926535*x[0]/L))',degree=1,L=lambd)

Xmax=lambd/2.
ddx = 0.01

# mesh
mesh = IntervalMesh(int(Xmax/ddx), 0, Xmax)

# basis functions
element = FiniteElement('Lagrange', interval, 1) #Lagrange polynomial of degree 1?

# mixed elements and associated function space
mixed_elements = MixedElement([element, element])
V = FunctionSpace(mesh, mixed_elements)

# Variables
hp = Function(V)
h, p = split(hp)
h_old = interpolate(h0, V.sub(0).collapse())

# Trial functions
u, v = split(TestFunction(V))

# Weak form of the equation to solve
eq1 = p*v*dx + inner(h.dx(0), v.dx(0))*dx
eq2 = h*u*dx - h_old*u*dx -dt*h**3*inner(p.dx(0),u.dx(0))*dx-dt/h*inner(h.dx(0),u.dx(0))*dx
eq = eq1+eq2

t=0
while t<t_max:
    t+=DT
    J = derivative(eq,hp)
    problem = NonlinearVariationalProblem(eq,hp,J=J)
    solver = NonlinearVariationalSolver(problem)
    prm = solver.parameters
    #prm["newton_solver"]["absolute_tolerance"] = 1E-6
    #prm["newton_solver"]["relative_tolerance"] = 1E-6
    #prm["newton_solver"]["maximum_iterations"] = 10000
    #prm["newton_solver"]["relaxation_parameter"] = 0.01
    
    # RANDOM INITIAL GUESS
    #hp.vector().set_local(rand(hp.vector().size()))
    #hp.vector().apply  
    
    solver.solve()
    
    height, height2 = hp.split(deepcopy=True)
    h_old.assign(height)
```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 26, 2020, 7:41pm UTC](https://fenicsproject.discourse.group/t/newton-solver-does-not-converge-due-to-non-linear-term/4176/2 "2020-09-26T19:41:16Z")

</div>

In your problem, if you do not start with a non-zero `hp`, the right hand side of your newton problem (Where F=eq1+eq2)  
`dF/d(hp)\Delta hp=-F(hp,uv)`  
Will be infinity, since eq(hp)=infinity.

I do not currently have the time to run your example. If I find time for it, i will keep you posted.

---

<div class="post-metadata">

**Author:** ![spoulain](https://avatars.discourse-cdn.com/v4/letter/s/0ea827/32.png) [@spoulain](https://fenicsproject.discourse.group/u/spoulain)\
**Post date:** [September 26, 2020, 8:18pm UTC](https://fenicsproject.discourse.group/t/newton-solver-does-not-converge-due-to-non-linear-term/4176/3 "2020-09-26T20:18:25Z")

</div>

I just searched more info taking into account your comment and apparently the default initial guess of the Newton’s solver is zero. I somehow overlooked that fact before so thanks, this is already quite helpful.

So I’m assuming i should either find a trick to remove the singularity in the equation, come up with a smarter initial guess, or use something else than Newton’s method. I’ll look into it.

edit: just taking another initial guess seems to be acceptable already, it was quite stupid… Good learning experience I guess. I’ll still need to deal with the singularity but that’s another issue