Linear Solvers

First we will consider the Poisson problem with $\partial_x u = 0$ on $\partial \Omega$. We do this by using $(1+\pi^2) \cdot cos(\pi\cdot x)$ as the rhs.

In [1]:
import import_hack
from methodsnm.mesh_2d import *
from methodsnm.fes import *

mesh = StructuredRectangleMesh(50, 50)
fes = P2_Triangle_Space(mesh)
source module for methodsNM imported.

As we have seen some times already, we construct the (Bi-)Linearforms with the according functions to solve the problem.

In [2]:
from methodsnm.forms import BilinearForm, LinearForm
from methodsnm.formint import *
from methodsnm.visualize import DrawFunction2D
from methodsnm.meshfct import *
from methodsnm.solver import *
from numpy import pi, cos

f = GlobalFunction(lambda x: (1+pi**2)*cos(pi*x[0]), mesh)
c = GlobalFunction(lambda x: 1, mesh=mesh)

blf = BilinearForm(fes)
blf += LaplaceIntegral(c)
blf += MassIntegral(c)
blf.assemble()
lf = LinearForm(fes)
lf += SourceIntegral(f)
lf.assemble()

uh = FEFunction(fes)
uh.vector = cg_solve(blf.matrix, lf.vector, n=250)
Warning: Using numerical differentiation for deriv evaluation in <class 'methodsnm.fe_2d.P2_Triangle_FE'> object.
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
In [3]:
DrawFunction2D(uh)
No description has been provided for this image

Comparison with Scipy Preconditioner

In [4]:
%timeit jacobi_prec(blf.matrix, len(lf.vector))
321 ms ± 969 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [5]:
from timeit import timeit
from scipy.sparse.linalg import gmres, cg

n = 250
%timeit gmres(blf.matrix, lf.vector, tol=1e-6, maxiter=n, M=jacobi_prec(blf.matrix, len(lf.vector)))
%timeit gmres(blf.matrix, lf.vector, tol=1e-6, maxiter=n)
%timeit cg(blf.matrix, lf.vector, tol=1e-6, maxiter=n, M=jacobi_prec(blf.matrix, len(lf.vector)))
%timeit cg(blf.matrix, lf.vector, tol=1e-6, maxiter=n)
%timeit cg_solve(blf.matrix, lf.vector, n=n, prec=True)
%timeit cg_solve(blf.matrix, lf.vector, n=n, prec=False)
836 ms ± 37.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
1.14 s ± 123 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
895 ms ± 32.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
96.6 ms ± 4.75 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
1.05 s ± 2.67 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
Reached limit of numbers of loop-execeutions! (in methodsnm.solver.cg_solve)
48.5 ms ± 329 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [6]:
from scipy.sparse.linalg import norm, inv
A = blf.matrix
print("unkonditioniertes A:", norm(A)*norm(inv(np.array(A))))
CA = np.matmul(jacobi_prec(A), A)
print("konditioniertes A:", norm(CA)*norm(inv(CA)))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[6], line 3
      1 from scipy.sparse.linalg import norm, inv
      2 A = blf.matrix
----> 3 print("unkonditioniertes A:", norm(A)*norm(inv(np.array(A))))
      4 CA = np.matmul(jacobi_prec(A), A)
      5 print("konditioniertes A:", norm(CA)*norm(inv(CA)))

File /usr/local/lib/python3.11/site-packages/scipy/sparse/linalg/_matfuncs.py:71, in inv(A)
     69 # Check input
     70 if not (scipy.sparse.issparse(A) or is_pydata_spmatrix(A)):
---> 71     raise TypeError('Input must be a sparse matrix')
     73 # Use sparse direct solver to solve "AX = I" accurately
     74 I = _ident_like(A)

TypeError: Input must be a sparse matrix
In [ ]: