【Scipy优化使用教程】六、寻根，不动点问题及其加速算法

本文链接：https://blog.csdn.net/w_w_y/article/details/123138697

参考官网：Scipy.

寻根，不动点问题及其加速算法

不动点问题
单变量方程求根
多变量方程求根（method=‘lm’）
大规模求根问题
加速方法

不动点问题

与寻找一个函数的零点密切相关的问题是寻找一个函数的不动点。一个函数的不动点是指对该函数进行求值时返回的点： $g (x) = x$ 。很明显，这是 $f (x) = g (x) - x$ 的零点问题。fixed_point提供了一个简单的迭代方法，使用Aitkens序列加速来估计 $g (x)$ 的不动点。

from scipy import optimize
import numpy as np
def func(x, c1, c2):
   return np.sqrt(c1/(x+c2))
c1 = np.array([10,12.])
c2 = np.array([3, 5.])
ans = optimize.fixed_point(func, [1.2, 1.3], args=(c1,c2))
print(ans)

结果为：

[1.4920333  1.37228132]

单变量方程求根

例如：
在这里插入图片描述

import numpy as np
from scipy.optimize import root
def func(x):
    return x + 2 * np.cos(x)
sol = root(func, 0.3)
print(sol.x)

其结果为：

[-1.02986653]

多变量方程求根（method=‘lm’）

例如：
在这里插入图片描述

import numpy as np
from scipy.optimize import root
def func2(x):
    f = [x[0] * np.cos(x[1]) - 4,
         x[1]*x[0] - x[1] - 5]
    df = np.array([[np.cos(x[1]), -x[0] * np.sin(x[1])],
                   [x[1], x[0] - 1]])
    return f, df
sol = root(func2, [1, 1], jac=True, method='lm')
print(sol.x)

其结果为：

[6.50409711 0.90841421]

大规模求根问题

hybr和 lm不能求解大规模问题，因为其要算N×N的雅可比矩阵，随着变量的增多，计算量会越来越大。
例如，考虑以下问题：我们需要解决以下微分方程：
在这里插入图片描述
限制条件为：上边缘满足 $P (x, 1) = 1$ ，边界上满足 $P = 0$ 。这可以通过用网格上的值来逼近连续函数P来实现，(个人理解就是流体用的有限差分法） $P_{n,m}≈P(nh,mh)$ 。然后可以对导数和积分进行近似处理；例如：
在这里插入图片描述
这样就可以使用大型求解器去求解：比如krylov, broyden2和 anderson。这些方法使用的是所谓的非精确牛顿法，它不是精确计算雅各布矩阵，而是对其形成一个近似值。

import numpy as np
from scipy.optimize import root
from numpy import cosh, zeros_like, mgrid, zeros

# 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 residual(P):
   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

# solve
guess = zeros((nx, ny), float)
sol = root(residual, guess, method='krylov', options={'disp': True})
#sol = root(residual, guess, method='broyden2', options={'disp': True, 'max_rank': 50})
#sol = root(residual, guess, method='anderson', options={'disp': True, 'M': 10})
print('Residual: %g' % abs(residual(sol.x)).max())

# visualize
import matplotlib.pyplot as plt
x, y = mgrid[0:1:(nx*1j), 0:1:(ny*1j)]
plt.pcolormesh(x, y, sol.x, shading='gouraud')
plt.colorbar()
plt.show()

结果为：

0:  |F(x)| = 40.1231; step 1
1:  |F(x)| = 16.821; step 1
2:  |F(x)| = 5.69025; step 1
3:  |F(x)| = 1.41135; step 1
4:  |F(x)| = 0.0386739; step 1
5:  |F(x)| = 0.00171956; step 1
6:  |F(x)| = 0.000132096; step 1
7:  |F(x)| = 6.2608e-06; step 1
8:  |F(x)| = 3.69188e-07; step 1
Residual: 3.69188e-07

在这里插入图片描述

加速方法

对于krylov方法来说，其大部分时间都用于计算雅可比矩阵的逆上面了。如果你有一个逆矩阵的近似值 $M≈J^{-1}$ ,那么你可以用它来对线性反演问题进行预处理。
我们不去求解 $J s = y$ 的问题，而是求解 $M J s = M y$ 。因为矩阵比原来更 "接近 "单位矩阵，所以Krylov方法应该更容易处理这个方程。
矩阵M可以作为一个选项options['jac_options']['inner_M']传递给Krylov方法。
对于上一节的问题，我们注意到要解决的函数由两部分组成：第一部分是拉普拉斯算子,第二个是积分。我们实际上可以很容易地计算出对应于拉普拉斯算子部分的雅各布系数：
在这里插入图片描述

$J_1$ 可以用 scipy.sparse.linalg.splu(或者用scipy.sparse.linalg.spilu近似)求出来。所以我们使用 $M≈J^{-1}_1$ 进行预处理，达到加速计算的效果。

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()