# Solving AU = F using linalg.cg results in 0 iterations

I am working on solving the following PDE: $$\left(\mu_{x}\frac{\partial^{2}u}{\partial x^{2}}+\mu_{y}\frac{\partial^{2}u}{\partial y^{2}}\right)=f(x,y) \tag 1$$

Which is then discretised:

$$- \mu_{x} \left(\frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{h^2}\right) - \mu_{y} \left(\frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{h^2}\right) = f_{i,j}$$

once all the nonlinear functions have been computed at the current guess $$u$$, we have $$A$$ as the Jacobian matrix in newtons method and $$F$$ for the right-hand side.

We now have a linear equation system $$JU = F$$. I am now trying to compute $$U$$ using sp.sparse.linalg.cg however, when I run it, it seems to give me 0 iterations and a result of 0. 0 usually indicates convergence but I can't see how it's managed to do this in 0 iterations. Would appreciate any insight into where I might be doing wrong with this.

I've provided a full reproducible example: solver.py

import scipy as sp
import numpy as np
from scipy.sparse.linalg import norm

def NewtonSys(fnon, jac, x0, tol, maxk, *fnonargs):
# fnon     - name of the nonlinear function f(x)
# jac      - name of the Jacobian function J(x)
# x0       - initial guess for the solution x0
# tol      - stopping tolerance for Newton's iteration
# maxk     - maximum number of Newton iterations before stopping
# fnonargs - optional arguments that will be passed to the nonlinear function (useful for additional function parameters)

k = 0
x = x0

F = eval(fnon)(x,*fnonargs)

Fsize = F.get_shape()
n = Fsize[0]

normF = (F.data ** 2).sum()

print(' k    f(xk)')

# Main Newton loop
while (normF > tol and k <= maxk):
# Evaluate Jacobian matrix
J = eval(jac)(x,n,fnon,F,*fnonargs)

# Take Newton step by solving the tangent problem
delta = np.linalg.solve(J.toarray(),F.toarray())
x = x - delta

F = eval(fnon)(x,*fnonargs)
normF = (F.data ** 2).sum()

print('{0:2.0f}  {1:2.2e}'.format(k, normF))
k += 1

if (k >= maxk):
print('Not converged')
else:
return x, J, F

def matrix( mu_x, mu_y, n ):

N = n**2

# We use a list-of-lists format for more efficient assembly
A = sp.sparse.lil_array((N, N), dtype=np.float64)
h = 1. / (n-1);

stencil = np.array([2 * mu_x + 2 * mu_y, mu_x * -1., mu_x * -1., mu_y * -1., mu_y * -1.]) / h**2

# Loop over each internal node in the grid, i,j = 1,2,...,n-2
# starts at the second row to avoid boundary and ends at n-2 to avoid boundary as well
for i in range(1,n-1):
for j in range(1,n-1):
# Find k-indices of the four neighbouring nodes
localStencilIndices = np.array([indexFD(i,j,n), indexFD(i+1,j,n), indexFD(i-1,j,n), indexFD(i,j+1,n), indexFD(i,j-1,n)])

# Add the local stencil for node (x_i,y_j) to the matrix
currentRow = indexFD(i,j,n)
for m in range (0,5):
A[currentRow, localStencilIndices[m]] = A[currentRow, localStencilIndices[m]] + stencil[m]

# After the matrix A has been assembled, we convert it to the column-major format
# for more efficient computations
A = A.tocsc()

return (A)

def indexFD(i, j, n):
return ( i*n + j )

import numpy as np

def boundaryConditions(n):

extNodes = []

# Find nodes at the boundary of the square
for j in range(0,n):
extNodes.append( indexFD(j,0,n) )
extNodes.append( indexFD(j,n-1,n) )
extNodes.append( indexFD(0,j,n) )
extNodes.append( indexFD(n-1,j,n) )

extNodes = np.unique(extNodes)
intNodes = np.setdiff1d(np.arange(0,n**2), extNodes);

return intNodes, extNodes

import scipy as sp
import copy

def fdJacobian(x,n,fnon,F0,*fnonargs):
# We replace the dense NumPy array with a sparse SciPy matrix
# (list-of-lists format is used for assembly)

J = sp.sparse.lil_array((n, n), dtype=np.float64)
h = 10e-8

for k in range(0,n):
xb = copy.deepcopy(x)
xb[k,0] = xb[k,0] + h

F = eval(fnon)(xb,*fnonargs)

for i in range(0,n):
J[i,k] = (F[i,0] - F0[i,0]) / h

# Return the Jacobian in column-major sparse format
return J.tocsc()

import math
def source_function(x, y, h):
x = x * h
y = y * h
if x >= 0.1 and x <= 0.3 and y >= 0.1 and y <= 0.3:
return 1
else:
return 0

def solve_rhs( u, A, n ):

N = n**2
capitalF = sp.sparse.lil_array((N, 1), dtype=np.float64)

# Find the internal nodes and the boundary nodes
intNodes, extNodes  = boundaryConditions(n)

Au = A[intNodes,:] @ u # Sparse matrix-vector multiplication
littleF = sp.sparse.lil_array((N, 1), dtype=np.float64)

for k in intNodes:
i, j = index_to_coordinates(k,n)
f_ij = source_function(i, j, 0.1)
littleF[k] = f_ij

capitalF[intNodes] = Au - littleF[intNodes]

# Set value of U at boundary nodes (this imposes u=0 at the boundary)
capitalF[extNodes] = u[extNodes]

return capitalF.tocsc()

def Solver(mu_x, mu_y, h):

n = int(1./h + 1)  # dimension of spatial mesh in each dimension
N = n**2      # total number of grid points

# Initial guess for Newton is the zero vector
u0 = sp.sparse.csc_array( (N,1), dtype=np.float64 )

# Assemble the Laplacian matrix
A = matrix(mu_x, mu_y, n)

# We plot the sparsity pattern of the matrix A to check that it has been
# correctly assembled.
fig = plt.figure()
cax = ax.matshow(A.toarray(), vmin=-4./h*mu_x, vmax=4./h*mu_y, cmap='coolwarm')
plt.title("Sparsity pattern")
fig.colorbar(cax)
plt.show()

# Newton iteration to solve nonlinear PDE
u, J, F = NewtonSys("solve_rhs", "fdJacobian", u0, 1e-6, 100, A, n )

return u, J, F

h = 0.1
#mu_x = 1.0, mu_y = 1.0
sol, J, F = Solver(1.0, 1.0, h)

# Solve AU = F
tolerance = 1e-6
max_iterations = 1000

def callback(xk):
callback.num_iterations += 1

callback.num_iterations = 0
x1, info1 = sp.sparse.linalg.cg(J, F.todense(), tol=tolerance, maxiter=max_iterations, callback=callback)
print("Iterations: ", callback.num_iterations)
print("Info: ", info1)

• why there is a newton step solving nonlinear PDE (see your comment)? Isn't just getting the matrix A and solve by cg? Dec 11, 2023 at 22:58
• There is a lot going on here, but it looks like you have both a finite difference function and a sparse array to compute the Laplacian being passed to a Newton solver for a linear problem. Try to reduce this to a minimal working example and I suspect you will figure out your problems. Dec 12, 2023 at 17:04