diff options
author | Ralf Gommers <ralf.gommers@gmail.com> | 2019-04-27 18:06:38 +0200 |
---|---|---|
committer | Ralf Gommers <ralf.gommers@gmail.com> | 2019-04-27 18:10:04 +0200 |
commit | cf0d554d90a4169a038778461122e8b07310d9e1 (patch) | |
tree | 6096dd841d23a6f0b76b8b5a31998da2c8500480 | |
parent | 6e1b9ff8580bb4aa1f49768301ba343d6b50ed36 (diff) | |
download | scipy-sphinx-theme-cf0d554d90a4169a038778461122e8b07310d9e1.tar.gz |
MAINT: fix most errors and warnings for `make html`
All that's left is some autodoc related warnings for scipy.odr
(from test_autodoc3.py)
-rw-r--r-- | examples/newton_krylov_preconditioning.py | 94 | ||||
-rw-r--r-- | test_optimize.rst | 5 |
2 files changed, 97 insertions, 2 deletions
diff --git a/examples/newton_krylov_preconditioning.py b/examples/newton_krylov_preconditioning.py new file mode 100644 index 0000000..4736157 --- /dev/null +++ b/examples/newton_krylov_preconditioning.py @@ -0,0 +1,94 @@ +import numpy as np +from scipy.optimize import root +from scipy.sparse import spdiags, kron +from scipy.sparse.linalg import spilu, LinearOperator +from numpy import cosh, zeros_like, mgrid, zeros, eye + +# parameters +nx, ny = 75, 75 +hx, hy = 1./(nx-1), 1./(ny-1) + +P_left, P_right = 0, 0 +P_top, P_bottom = 1, 0 + +def get_preconditioner(): + """Compute the preconditioner M""" + diags_x = zeros((3, nx)) + diags_x[0,:] = 1/hx/hx + diags_x[1,:] = -2/hx/hx + diags_x[2,:] = 1/hx/hx + Lx = spdiags(diags_x, [-1,0,1], nx, nx) + + diags_y = zeros((3, ny)) + diags_y[0,:] = 1/hy/hy + diags_y[1,:] = -2/hy/hy + diags_y[2,:] = 1/hy/hy + Ly = spdiags(diags_y, [-1,0,1], ny, ny) + + J1 = kron(Lx, eye(ny)) + kron(eye(nx), Ly) + + # Now we have the matrix `J_1`. We need to find its inverse `M` -- + # however, since an approximate inverse is enough, we can use + # the *incomplete LU* decomposition + + J1_ilu = spilu(J1) + + # This returns an object with a method .solve() that evaluates + # the corresponding matrix-vector product. We need to wrap it into + # a LinearOperator before it can be passed to the Krylov methods: + + M = LinearOperator(shape=(nx*ny, nx*ny), matvec=J1_ilu.solve) + return M + +def solve(preconditioning=True): + """Compute the solution""" + count = [0] + + def residual(P): + count[0] += 1 + + d2x = zeros_like(P) + d2y = zeros_like(P) + + d2x[1:-1] = (P[2:] - 2*P[1:-1] + P[:-2])/hx/hx + d2x[0] = (P[1] - 2*P[0] + P_left)/hx/hx + d2x[-1] = (P_right - 2*P[-1] + P[-2])/hx/hx + + d2y[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2])/hy/hy + d2y[:,0] = (P[:,1] - 2*P[:,0] + P_bottom)/hy/hy + d2y[:,-1] = (P_top - 2*P[:,-1] + P[:,-2])/hy/hy + + return d2x + d2y + 5*cosh(P).mean()**2 + + # preconditioner + if preconditioning: + M = get_preconditioner() + else: + M = None + + # solve + guess = zeros((nx, ny), float) + + sol = root(residual, guess, method='krylov', + options={'disp': True, + 'jac_options': {'inner_M': M}}) + print('Residual', abs(residual(sol.x)).max()) + print('Evaluations', count[0]) + + return sol.x + +def main(): + sol = solve(preconditioning=True) + + # visualize + import matplotlib.pyplot as plt + x, y = mgrid[0:1:(nx*1j), 0:1:(ny*1j)] + plt.clf() + plt.pcolor(x, y, sol) + plt.clim(0, 1) + plt.colorbar() + plt.show() + + +if __name__ == "__main__": + main() diff --git a/test_optimize.rst b/test_optimize.rst index a09e0d7..0a55c28 100644 --- a/test_optimize.rst +++ b/test_optimize.rst @@ -29,7 +29,7 @@ The module contains: 5. Multivariate equation system solvers (:func:`root`) using a variety of algorithms (e.g. hybrid Powell, Levenberg-Marquardt or large-scale - methods such as Newton-Krylov). + methods such as Newton-Krylov [KK]_). Below, several examples demonstrate their basic usage. @@ -771,7 +771,8 @@ lot more depth to this topic than is shown here. .. rubric:: References -Some further reading and related software: +Some further reading and related software for solving large-scale problems +( [PP]_, [AMG]_): .. [KK] D.A. Knoll and D.E. Keyes, "Jacobian-free Newton-Krylov methods", J. Comp. Phys. 193, 357 (2003). |