I have a nonlinear function fun with n inputs and n outputs. I also have a function jac which calculates the Jacobian, which is sparse. Which SciPy method, if any, can I use to find the root of fun?

I am surprised that I can't find a method in scipy.optimize that takes a Jacobian in sparse matrix format. But I believe this can be simply solved by the Newton's method, and so I am confused if such a method does not exist. There are many quasi-Newton methods in this package that estimate the Jacobian, but they do not seem quite right. The only thing that looks similar is newton_krylov which documents something about sparse preconditioners. Is this what I should be looking at? If there is no such method, should I just do something like the following?

while not some_stopping_condition(x):
    x -= scipy.sparse.spsolve(jac(x), fun(x))

I would prefer to use an existing SciPy function in case I do something stupid. Notably, as it happens, I believe my particular problem is poorly conditioned, and I think I would benefit from being able to use scipy.sparse.lsqr or lsmr which takes a damping parameter.

Incidentally, the problem is solving large-displacement elasticity equations on a 2D membrane in 3D space with FEM. The problem is sparse because it is FEM. The problem is nonlinear because of the large displacements (and the corresponding Green-Lagrange strains). The poor conditioning comes from the fact that a small perturbation of any node in the direction perpendicular to the surface creates little change in strain energy (and the flatter the surface, the smaller the change).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.