Calculating condition number in a nonlinear problem

variational formulation
1

Hello
I want to calculate the condition number of stiffness matrix in nonlinear Poisson equation. I followed this post. Here is the minimal example:

from dolfin import *
import numpy as np

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1)

# Define boundary condition
g = Constant(1.0)
bc = DirichletBC(V, g, DirichletBoundary())

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])",degree=1)
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx

# Compute solution
solve(F == 0, u, bc)

A = assemble(lhs(F))
print (np.linalg.cond(A.array()))

I am expecting to get the condition number of the stiffness matrix but it returns this error:

    print (np.linalg.cond(A.array()))
AttributeError: 'float' object has no attribute 'array'

I am just wondering how I can get it done. Any help is appreciated.
Best

2

When you have a non linear problem, you use Newtons method to create linear systems to solve (as F is a vector).
This means that you are solving the problem

J=-F(u,v)

where J=derivative(F, u, TrialFunction(V)).
Note that the matrix J is dependent on the current solution of u i.e. the condition number changes during the solve process.
See for instance: https://fenicsproject.org/qa/2913/newtons-method-programmed-manually-for-parabolic-pde/

1 Like
3

In addition to @dokken’s comment, you can use a NewtonSolver and corresponding implementation of a NonlinearProblem to get the linearised operator at each Newton iteration. This is probably a good spot to compute (or estimate) the condition number. See, e.g. here.

To clarify: in each Newton iteration we compute updates to the estimate of our solution given some known initial guess u^{0}_h

u^{m+1}_h = u^{m}_h + \delta u^m_h \quad m = 0, 1, 2, \ldots.

The update \delta u^m_h is the solution of the following system: find \delta u^m_h \in V^h such that

J_u (\delta u^{m}_h; u^{m}_h, v) = - F(u^{m}_h, v) \quad \forall v \in V^h

where F(u, v) is the residual formulation, J_u = F^\prime(u, v) is its Gateaux derivative and V^h is some appropriate finite element space.

Note here it’s the form J_u(\cdot; \cdot, \cdot) which is bilinear (\delta u^m_h is unknown and u^m_h is known) and therefore its discretisation by finite elements yields the underlying linear operator. To write lhs(F) as in your code does not make sense since the residual is a linear form (and therefore a vector).

2 Likes
4

Thanks for your clarifications. I learned that condition number changes during the solve process. With that being said, what should I add to the above example code to calculate the condition number in each iteration (or the last iteration)? I just want to see it in a practical code.