# Newton method programmed manually

Hi all, This is probably a simple thing but I've been hitting my head against it. I need to program the newton method manually, and want to compare it to the nonlinearvariationalsolver. Both converge but I can't seem to match my residual norms. Here is my code:

from dolfin import *
import numpy as np

mesh = IntervalMesh(100, 0,1)
V = FunctionSpace(mesh, "Lagrange",1)    # Order 1 function space

dU = TrialFunction(V)
test_U = TestFunction(V)
U = Function(V)

left, right = compile_subdomains([
"(std::abs( x[0] )           < DOLFIN_EPS) && on_boundary",
"(std::abs( x[0]-1.0 ) < DOLFIN_EPS) && on_boundary"])

bcs = [
DirichletBC(V, 1, left),
DirichletBC(V, 0 , right),
]

k = 1.0*U+1.0

Jac = derivative(F, U, dU)

U_init = Expression("1-x[0]")

U.interpolate(U_init)

if False:
problem = NonlinearVariationalProblem(F, U, bcs, Jac) #,form_compiler_parameters=ffc_options
solver = NonlinearVariationalSolver(problem)
solver.solve()
else:
a_tol, r_tol = 1e-7, 1e-10
bcs_du = homogenize(bcs)
U_inc = Function(V)
nIter = 0
eps = 1

while eps > 1e-10 and nIter < 10:              # Newton iterations
nIter += 1
A, b = assemble_system(Jac, -F, bcs_du)
solve(A, U_inc.vector(), b)     # Determine step direction
eps = np.linalg.norm(U_inc.vector().array(), ord = 2)

a = assemble(F)
fnorm = np.linalg.norm(a.array(), ord = 2)
lmbda = 1.0     # step size, initially 1

U.vector()[:] += lmbda*U_inc.vector()    # New u vector

print '      {0:2d}  {1:3.2E}  {2:5e}'.format(nIter, eps, fnorm)

plot(U)
interactive()


Change the True / False above to switch methods. With the nonlinearvariationalsolver I get the output:

Solving nonlinear variational problem.
Newton iteration 1: r (abs) = 5.898e-03 (tol = 1.000e-10) r (rel) = 5.927e-02 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 6.381e-06 (tol = 1.000e-10) r (rel) = 6.413e-05 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 7.288e-12 (tol = 1.000e-10) r (rel) = 7.325e-11 (tol = 1.000e-09)
Newton solver finished in 3 iterations and 3 linear solver iterations.


With my code I get:

   1  6.24E-01  2.236057e+00
2  1.56E-02  2.121959e+00
3  1.14E-05  2.121320e+00
4  7.20E-12  2.121320e+00


Can anyone point out my error?

Thanks,
Mike

In your second column there is a norm of the increment while NewtonSolver checks convergence using norm of the equation residual by default. You can change this behaviour by

solver.parameters['newton_solver']['convergence_criterion'] = 'incremental'


Thanks Jan, but I am interested in how to calculate the norm of the equation residual actually. I thought thats what I was doing with the commands:

a = assemble(F)
fnorm = np.linalg.norm(a.array(), ord = 2)


But it doesn't seem to match the output from the newton solver... Maybe I have teh wronge equation?

There has been a little update on how and when the norms are calculated in Dolfin's Newton method, cf. https://bitbucket.org/fenics-project/dolfin/commits/0345014b7cd09bff6bbd8f56665a2832b13ddebb. You may want to build from Git master.

Nico, this is barely significant to the problem. The problem actually is that mwelland's home-made Newton solver does not show convergence in terms of $||F||_{\ell^2}$. He gets still $||F||_{\ell^2}\approx 2$

Thanks for the comments guys, looks like my problem was not applying the bcs_du to the matrix 'a' before calculating the norm. I tried just b.norm (inspired from nschole's notes) and it matches the nonlinearvariationalsolver output.

answered Jul 14, 2013 by FEniCS User (8,410 points)
selected Jul 14, 2013

Yeah, it was tricky for me too to realize that assemble(F) needs application of bcs_du because test space of Dirichlet problem F is constrained by bcs_du ;-)

can you clarify exactly what you did to make it work, i.e. code?

I guess

...

a = assemble(F)
fnorm = norm(a, 'l2') # better than numpy function
lmbda = 1.0     # step size, initially 1

...


actually, I just replaced the calculation of fnorm to
fnorm = b.norm('l2')
(which avoids having to assemble F again too!) Here is modified code with a bit more description:

from dolfin import *
import numpy as np

mesh = IntervalMesh(100, 0,1)
V = FunctionSpace(mesh, "Lagrange",1)    # Order 1 function space

dU = TrialFunction(V)
test_U = TestFunction(V)
U = Function(V)

left, right = compile_subdomains([
"(std::abs( x[0] )           < DOLFIN_EPS) && on_boundary",
"(std::abs( x[0]-1.0 ) < DOLFIN_EPS) && on_boundary"])

bcs = [
DirichletBC(V, 1, left),
DirichletBC(V, 0 , right),
]

k = 1.0*U+1.0

Jac = derivative(F, U, dU)

U_init = Expression("1-x[0]")

U.interpolate(U_init)

if False:
problem = NonlinearVariationalProblem(F, U, bcs, Jac) #,form_compiler_parameters=ffc_options
solver = NonlinearVariationalSolver(problem)
solver.solve()
else:
a_tol, r_tol = 1e-7, 1e-10
bcs_du = homogenize(bcs)
U_inc = Function(V)
nIter = 0
eps = 1

while eps > 1e-10 and nIter < 10:              # Newton iterations
nIter += 1
A, b = assemble_system(Jac, -F, bcs_du)
solve(A, U_inc.vector(), b)     # Determine step direction
eps = np.linalg.norm(U_inc.vector().array(), ord = 2)

a = assemble(F)
for bc in bcs_du:
bc.apply(a)
print 'b.norm("l2") = ', b.norm('l2'), 'np.linalg.norm(a.array(), ord = 2) = ', np.linalg.norm(a.array(), ord = 2)
fnorm = b.norm('l2')
lmbda = 1.0     # step size, initially 1

U.vector()[:] += lmbda*U_inc.vector()    # New u vector

print '      {0:2d}  {1:3.2E}  {2:5e}'.format(nIter, eps, fnorm)

plot(U)
interactive()


Hope this helps!

perfect, thanks!

I'm having an issue with the solver when solving a system of equations :

from dolfin import *
from pylab  import plot, show, spy
import sys

bb       = bool(int(sys.argv[1]))

mesh     = IntervalMesh(4, 0, 1)
Q        = FunctionSpace(mesh, "Lagrange",1)       # Order 1 function space
MQ       = Q * Q

dU       = TrialFunction(MQ)
dU1, dU1 = split(dU)
test_U   = TestFunction(MQ)
t1, t2   = split(test_U)
U        = Function(MQ)
u1, u2   = split(U)

left, right = compile_subdomains([
"(std::abs( x[0] )     < DOLFIN_EPS) && on_boundary",
"(std::abs( x[0]-1.0 ) < DOLFIN_EPS) && on_boundary"])

bcs = [DirichletBC(MQ.sub(0), 1, left),
DirichletBC(MQ.sub(0), 0, right),
DirichletBC(MQ.sub(1), 0, left),
DirichletBC(MQ.sub(1), 1, right)]

k   = 1.0*u1 + 1.0

vnorm    = sqrt(dot(u2, u2) + 1e-10)
cellh    = CellSize(mesh)
t2hat    = t2 + cellh/(2*vnorm)*dot(u2, t2.dx(0))

F1  = + u1 * u1.dx(0) * t1.dx(0) * dx \
+ u1.dx(0) * t1.dx(0) * dx
c   = 1e-9
F2  = c * u2 * t2hat * dx
#F2  = u2 * t2 * dx
F   = F1 + F2

u1_init = Expression("1-x[0]")
u2_init = Constant(1)

U_k = project(as_vector([u1_init, u2_init]), MQ) # project inital values
U.vector().set_local(U_k.vector().array())       # initalize u1, u2 in solution

Jac     = derivative(F, U, dU)

if bb:
problem = NonlinearVariationalProblem(F, U, bcs, Jac)
solver  = NonlinearVariationalSolver(problem)
solver.solve()
else:
atol, rtol = 1e-7, 1e-10                      # abs/rel tolerances
lmbda      = 1.0                              # relaxation parameter
U_inc      = Function(MQ)                     # residual
bcs_u      = homogenize(bcs)                  # residual is zero on boundary
nIter      = 0                                # number of iterations
residual   = 1                                # residual
rel_res    = residual                         # relative residual
maxIter    = 100                              # max iterations

while residual > atol and rel_res > rtol and nIter < maxIter:
nIter  += 1                                 # increment interation
A, b    = assemble_system(Jac, -F, bcs_u)   # assemble system
solve(A, U_inc.vector(), b)                 # determine step direction
rel_res = U_inc.vector().norm('l2')         # calculate norm

# calculate absolute residual :
a = assemble(F)
for bc in bcs_u:
bc.apply(a)
residual = b.norm('l2')

U.vector()[:] += lmbda*U_inc.vector()       # New u vector

string = "Newton iteration %d: r (abs) = %.3e (tol = %.3e) " \
+"r (rel) = %.3e (tol = %.3e)"
print string % (nIter, residual, atol, rel_res, rtol)

xcrd = mesh.coordinates()[:,0]
plot(xcrd, project(u1, Q).vector().array())
plot(xcrd, project(u2, Q).vector().array())

show()


it appears that the u2 solution's boundary conditions are not being applied properly. Any ideas?