python求解偏微分方程

python求解偏微分方程

偏微分方程（Partial Differential Equations, PDE）是研究连续介质中的许多物理现象所必需的重要数学工具。PDE 涉及了空间、时间和其它同步变量之间的关系，因此对于有限元分析（FEM）和流体力学等领域来说，具有极为重要的应用价值。下面我们将简单介绍使用 Python 求解偏微分方程的基本方法。

1. 引入库

在 Python 中，我们可以使用 SciPy 和 NumPy 库来处理偏微分方程。其中，NumPy 用于数值计算，而 SciPy 则提供了一些特定的算法，包括线性方程组求解、优化、数值积分和微分方程等。因此，我们需要在程序中引入这两个库：

```python

import numpy as np

from scipy import sparse

from scipy.sparse.linalg import spsolve

```

2. 构建矩阵

在求解偏微分方程时，我们通常需要构建雅可比矩阵。这里举一个简单的例子，设有一个一维热传导方程：

$$ \frac{\partial^2 u}{\partial x^2} = f(x) $$

其中，$u$ 是未知函数，$f(x)$ 是给定函数。为了求解这个方程，我们可以按照离散化的方法来处理。

我们将区间 $[0,1]$ 分成 $n$ 个小区间，即 $x_0 = 0$，$x_n

= 1$，$x_i = ih$，$h = 1/n$。因此，$u(x_i)$ 可以用 $u_i$ 来表示。我们将前式中的二阶导数离散化，得到如下近似式：

$$ \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} = f_i $$

解出 $u_i$，得到：

$$u_i = \frac{1}{h^2} \left(u_{i+1} + u_{i-1} - h^2

f_i\right)$$

这样，我们就可以得到一个线性方程组：

$$A\mathbf{u} = \mathbf{f}$$

其中，$\mathbf{u}$ 是 $[u_1, u_2, ..., u_n]$ 的列向量，

$\mathbf{f}$ 是 $[f_1, f_2, ..., f_n]$ 的列向量。

对于边界条件，我们通常需要指定 $u_0$ 和 $u_n$ 的值。这里

我们以 $u_0 = 0$，$u_n = 1$ 为例。

下面的代码演示了如何构建这个矩阵：

```python

n = 50

h = 1/n

x = np.linspace(0, 1, n+1)

A = np.zeros((n+1, n+1))

for i in range(1, n):

A[i,i-1] = 1

A[i,i] = -2

A[i,i+1] = 1

A[0,0] = 1

A[n,n] = 1

f = (x)

f[0] = 0

f[n] = 1

```

这里我们使用了 NumPy 的 linspace 函数生成了一个长度为

$n+1$ 的等差数组，用来存储 $x_i$。然后，我们创建了一个

$n+1$ 行 $n+1$ 列的零矩阵 $A$，使用 for 循环填充了矩阵的非对角线部分。最后，我们给矩阵的第一行和最后一行赋值，以满足边界条件。同时，我们还生成了一个长度为 $n+1$ 的向量 $f$，用于存储$[f_1, f_2, ..., f_n]$。

3. 求解线性方程组

得到线性方程组之后，我们还需要求解它的解。这里我们使用Scipy 提供的函数 spsolve 来求解。代码如下：

```python

u = spsolve(sparse.csr_matrix(A), f)

```

在这里，我们首先将 NumPy 数组转换成 Scipy 稀疏矩阵，然后使用 spsolve 函数求解线性方程组。得到的解 $u$ 是一个长度为$n+1$ 的向量，其中 $u_i$ 表示 $u(x_i)$ 的近似值。

4. 可视化

最后，我们可以使用 Matplotlib 库来可视化求解结果：

```python

import matplotlib.pyplot as plt

(x, u, 'r-', label='Numerical')

(x, (x), 'g--', label='Exact')

plt.xlabel('x')

plt.ylabel('u(x)')

(True)

plt.legend()

()

```

这里我们将近似解 $u(x)$ 画成红色线条，将准确解

$\sin(x)$ 画成绿色虚线。这样我们就可以直观地看到我们的数值解与准确解之间的关系。

完整代码如下：

```python

import numpy as np

from scipy import sparse

from scipy.sparse.linalg import spsolve

import matplotlib.pyplot as plt

n = 50

h = 1/n

x = np.linspace(0, 1, n+1)

A = np.zeros((n+1, n+1))

for i in range(1, n):

A[i,i-1] = 1