在 scipy 中求解离散边界值问题

日期2016-10-22（最后修改）

数学公式¶

我们考虑一个在正方形域 $\Omega = [0, 1] \times [0, 1]$ 上的非线性椭圆边界值问题：$$ \Delta u + k f(u) = 0 \\ u = 0 \text{ on } \partial \Omega $$

这里 $u = u(x,y)$ 是一个未知函数，$\Delta$ 是 拉普拉斯算子，$k$ 是某个常数，$f(u)$ 是一个给定函数。

对这类问题，通常的计算方法是离散化。我们使用一个步长为 $h$ 的均匀网格，并定义 $u_{i, j} = u(i h , jh)$。通过使用 5 点有限差分 来近似拉普拉斯算子，我们得到一个方程组，形式如下：$$ u_{i - 1, j} + u_{i + 1, j} + u_{i, j - 1} + u_{i, j + 1} - 4 u_{i,j} + c f(u_{i,j}) = 0 \\ u_{i, j} = 0 \text{ on } \partial \Omega $$

这里 $c = k h^2$。

为 scipy 定义问题¶

从现在开始，我们关注离散版本，并考虑一个每个维度有 $n = 100$ 个刻度的网格，并设置 $c = 1$ 和 $f(u) = u^3$。

在 [1]

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import least_squares
from scipy.sparse import coo_matrix

在 [2]

n = 100
c = 1

在 [3]

def f(u):
    return u**3

def f_prime(u):
    return 3 * u**2

为了求解方程组，我们将使用 scipy.optimize.least_squares。

我们定义一个函数来计算每个方程的左侧。注意，我们假设边界上的值为零，并且在优化过程中不会改变它们。

在 [4]中

def fun(u, n, f, f_prime, c, **kwargs):
    v = np.zeros((n + 2, n + 2))
    u = u.reshape((n, n))
    v[1:-1, 1:-1] = u
    y = v[:-2, 1:-1] + v[2:, 1:-1] + v[1:-1, :-2] + v[1:-1, 2:] - 4 * u + c * f(u)
    return y.ravel()

始终建议尽可能提供解析雅可比矩阵。在我们的问题中，我们有 $n^2=10000$ 个方程和变量，但每个方程仅依赖于少数变量，因此我们应该以稀疏格式计算雅可比矩阵。

预先计算雅可比矩阵中非零元素的行和列索引很方便。我们定义了相应的函数

在 [5]中

def compute_jac_indices(n):
    i = np.arange(n)
    jj, ii = np.meshgrid(i, i)

    ii = ii.ravel()
    jj = jj.ravel()

    ij = np.arange(n**2)

    jac_rows = [ij]
    jac_cols = [ij]

    mask = ii > 0
    ij_mask = ij[mask]
    jac_rows.append(ij_mask)
    jac_cols.append(ij_mask - n)

    mask = ii < n - 1
    ij_mask = ij[mask]
    jac_rows.append(ij_mask)
    jac_cols.append(ij_mask + n)

    mask = jj > 0
    ij_mask = ij[mask]
    jac_rows.append(ij_mask)
    jac_cols.append(ij_mask - 1)

    mask = jj < n - 1
    ij_mask = ij[mask]
    jac_rows.append(ij_mask)
    jac_cols.append(ij_mask + 1)

    return np.hstack(jac_rows), np.hstack(jac_cols)

之后，在 coo_matrix 格式中计算雅可比矩阵很简单

在 [6]中

jac_rows, jac_cols = compute_jac_indices(n)

在 [7]中

def jac(u, n, f, f_prime, c, jac_rows=None, jac_cols=None):
    jac_values = np.ones_like(jac_cols, dtype=float)
    jac_values[:n**2] = -4 + c * f_prime(u)
    return coo_matrix((jac_values, (jac_rows, jac_cols)), shape=(n**2, n**2))

解决问题¶

在没有对问题的任何了解的情况下，我们最初将所有值设置为 0.5。请注意，不能保证给定的连续或离散问题具有唯一解。

在 [8]中

u0 = np.ones(n**2) * 0.5

预先计算雅可比矩阵中非零元素的行和列

在 [9]中

jac_rows, jac_cols = compute_jac_indices(n)

现在我们准备运行优化。第一个解将在不对方程施加任何约束的情况下计算。

在 [10]中

res_1 = least_squares(fun, u0, jac=jac, gtol=1e-3, args=(n, f, f_prime, c), kwargs={'jac_rows': jac_rows, 'jac_cols': jac_cols}, verbose=1)

`gtol` termination condition is satisfied.
Function evaluations 106, initial cost 1.0412e+02, final cost 5.2767e-03, first-order optimality 9.04e-04.

下面我们可视化第一个解。左侧图显示了扁平化的解，中间图显示了解在方形域中的外观，右侧图显示了每个节点的最终残差。

在 [11]中

plt.figure(figsize=(16, 5))
plt.subplot(132)
plt.imshow(res_1.x.reshape((n, n)), cmap='coolwarm', vmin=-max(abs(res_1.x)), vmax=max(abs(res_1.x)))
plt.colorbar(use_gridspec=True, fraction=0.046, pad=0.04)
plt.subplot(131)
plt.plot(res_1.x)
plt.subplot(133)
plt.plot(res_1.fun)
plt.tight_layout()

某些物理考虑可能要求解在任何地方都必须非负。我们可以通过向求解器指定约束来实现这一点

在 [12]中

res_2 = least_squares(fun, u0, jac=jac, bounds=(0, np.inf), gtol=1e-3,
                      args=(n, f, f_prime, c), kwargs={'jac_rows': jac_rows, 'jac_cols': jac_cols}, verbose=1)

`gtol` termination condition is satisfied.
Function evaluations 34, initial cost 1.0412e+02, final cost 4.1342e-02, first-order optimality 9.55e-04.

在 [13]中

plt.figure(figsize=(16, 5))
plt.subplot(132)
plt.imshow(res_2.x.reshape((n, n)), cmap='Reds')
plt.colorbar(use_gridspec=True, fraction=0.046, pad=0.04)
plt.subplot(131)
plt.plot(res_2.x)
plt.subplot(133)
plt.plot(res_2.fun)
plt.tight_layout()

我们看到设置下界允许我们找到不同的非负解。

您可以尝试从不同的起点运行优化，使用不同的约束或更改 $c$ 和 $f(u)$，并查看它如何影响 least_squares 求解器的结果。

在 [ ]

章节作者：Nikolay Mayorov