I’m currently working on solving the following PDE: $$\begin{equation} -(\mu_x \frac{\partial^2 u}{\partial x^2} + \mu_y \frac{\partial^2 u}{\partial y^2}) = f(x, y)\end{equation}$$

Where a right hand side F is formed, which we can use to solve Ax = b (AU = F) I’m using Python and the scipy.sparse.linalg.spilu function for the ILU decomposition to pass L and U to an iterative method to approximate the solution.

import scipy as sp
import numpy as np

# Jacobian from newton, will be left hand side A
u, J = NewtonSys("solve_f", "fdJacobian", u0, 1e-6, 100, A, n)
drop_tol = 0.1
LU = sp.sparse.linalg.spilu(J, drop_tol)
L = LU.L
U = LU.U

When I run this code, I encounter a singular matrix error during the ILU decomposition. This is unexpected because the Jacobian matrix J should not be singular.

The Jacobian matrix J is defined as follows: ${\bf u}_{int}$ are the values in the interior nodes and ${\bf u}_{bdary}$ the values in the boundary nodes. The Jacobian matrix J is defined as follows:

enter image description here

where ${\bf u}{int}$ are the values at the interior nodes, ${\bf u}{bdary}$ are the values at the boundary nodes, A is the diffusion matrix, and I is the identity matrix

I've tried adjusting the drop tolerance drop_tol in the spilu function and the h value in the finite difference approximation used to calculate J, but neither of these changes resolved the issue. Does anyone have any suggestions on what might be causing this error and how to fix it? Any help would be greatly appreciated!

Full solver for reproducible example as the parts need to work together for the problem to make sense. full_solver.py

  • $\begingroup$ You assume that scipy was wrong in telling you that the matrix is singular. But have you checked that it is indeed not singular? $\endgroup$ Dec 14, 2023 at 6:13


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.