计算科学 - 圆上的有限差分离散化 - 吾爱随笔录

圆上的有限差分离散化

计算科学有限差分离散化数字

2021-12-20 19:21:16

我正在尝试离散化微分算子 $\frac{d^2}{dx^2}$ 作用于 $S^1 = [0,1]$ 在一个圆周围使用有限多个点 $0, \frac{1}{N}, \frac{2}{N}, \dots, \frac{N-1}{N}$ .

这是矩阵 $N=5$ ：

[\begin{array}{cc} - 2 & 1 & 0 & 0 & 1 \\ 1 & - 2 & 1 & 0 & 0 \\ 0 & 1 & - 2 & 1 & 0 \\ 0 & 0 & 1 & - 2 & 1 \\ 1 & 0 & 0 & 1 & - 2 \end{array}]

$\left[ \begin{array}{cc} -2 & 1 & 0 & 0 & 1 \\ 1 & -2 & 1 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 1 & -2 & 1 \\\ 1 &0 & 0 & 1 & -2\end{array} \right]$

如果我们解方程 $\frac{d^2 f}{dx^2} = \lambda f$ 我们应该得到 $f(x) = e^{\lambda x}$ 和 $\lambda = 2\pi k$ 和 $k^2 \in \mathbb{Z}$ .

当我在 numpy 中实现这个时，重新缩放似乎都是错误的。有人可以解释一下这背后的理论吗？

N = 100
A = np.zeros((N,N))
d = np.arange(N)
A[d,d] = -2
A[d,(d+1)%N]=1
A[d,(d-1)%N]=1
A

生成的矩阵就像我说的那样。这里 $N=100$ .

array([[-2.,  1.,  0., ...,  0.,  0.,  1.],
       [ 1., -2.,  1., ...,  0.,  0.,  0.],
       [ 0.,  1., -2., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ..., -2.,  1.,  0.],
       [ 0.,  0.,  0., ...,  1., -2.,  1.],
       [ 1.,  0.,  0., ...,  0.,  1., -2.]])

我想出了这个针对特征值的临时缩放 $-A$ 将它们除以 $(2 \pi N)^2$ .

L = np.linalg.eig(-A)[0]/(2*np.pi)**2*N**2
L = np.round(L)
np.sort(L).astype(int)

结果看起来非常接近。我得到了完美正方形的序列。这是否意味着 $\frac{1}{2\pi N}A$ 是正确的重新调整？

array([   0,    1,    1,    4,    4,    9,    9,   16,   16,   25,   25,
         36,   36,   48,   48,   63,   63,   79,   79,   97,   97,  116,
        116,  137,  137,  160,  160,  184,  184,  209,  209,  235,  235,
        263,  263,  291,  291,  320,  320,  350,  350,  381,  381,  412,
        412,  443,  443,  475,  475,  507,  507,  538,  538,  570,  570,
        602,  602,  633,  633,  663,  663,  693,  693,  722,  722,  751,
        751,  778,  778,  804,  804,  830,  830,  853,  853,  876,  876,
        897,  897,  916,  916,  934,  934,  951,  951,  965,  965,  978,
        978,  988,  988,  997,  997, 1004, 1004, 1009, 1009, 1012, 1012,
       1013])

1个回答

您用于特征值和特征函数的表达式是错误的（根据 Wolfgang Bangerth 的评论）；因此，您得到的结果根本没有意义。

对于连续和离散情况，特征值都有解析表达式 $d^2/dx^2$ 在一个圆形段上——周期性的 BCs。

对于一个连续的情况 $S^1=[0,L]$ , $j$ 特征值 $\lambda_j$ ：

λ_{j}^{cont} = {\begin{cases} - \frac{π^{2}}{L^{2}} j^{2}, & j is even \\ - \frac{π^{2}}{L^{2}} (j - 1)^{2}, & j is odd \end{cases}

$\lambda_j^{\text{cont}} = \begin{cases} -\frac{\pi^2}{L^2}j^2,&\quad j \text{ is even}\\ -\frac{\pi^2}{L^2}(j-1)^2,&\quad j \text{ is odd} \end{cases}$

对于离散情况 $N$ 点均匀分布在一个圆和中心差分方案周围：

λ_{j}^{disc} = {\begin{cases} - \frac{4}{h^{2}} \sin^{2} (\frac{π j^{2}}{2 N}), & j is even \\ - \frac{4}{h^{2}} \sin^{2} (\frac{π (j - 1)^{2}}{2 N}), & j is odd \end{cases}

$\lambda_j^{\text{disc}} = \begin{cases} -\frac{4}{h^2}\sin^2\left(\frac{\pi j^2}{2N}\right),&\quad j \text{ is even}\\ -\frac{4}{h^2}\sin^2\left(\frac{\pi (j-1)^2}{2N}\right),&\quad j \text{ is odd} \end{cases}$

这实际上是矩阵特征值的解析表达式 $A$ . 这里， $h=L/N.$ 对于圆段 $S^1$ .

我编写了一个简单的 Python 脚本，用于比较连续运算符和离散运算符的特征值（使用来自实际矩阵numpy.linalg.eigvals）。

为了 $N=5$ ，第一个特征值的相对误差为：

[0.1248598 0.1248598 0.4272133 0.4272133]（第一个特征值为0并被丢弃，特征值的重数为2）

为了 $N=100$ ,

[0.00032894 0.00032894 0.00131525 0.00131525]

因此，该方案“收敛”。

#!/usr/local/bin/python
import numpy as np

# defining problem
L = 1.      # length of the circular interval S = [0,L] (starting is assumed to be fixed
N = 100       # Number of points discretizing circular interval S
h = L/N     # FD step-size (note, it is a circle; thus, N vs. N-1

# forming discretized d^2 / dx^2 operator
A = np.zeros((N, N))
d = np.arange(N)
A[d, d] = -2
A[d, (d+1) % N] = 1
A[d, (d-1) % N] = 1
A *= 1./(h**2)

# calculating eigenvalues
A_eig = np.linalg.eigvals(A)
A_eig = np.flip(np.sort(A_eig), 0)  # decreasing order

# comparing to analytical expressions for eigenvalues in continuous and discrete cases
eig_cont_anal = np.zeros(N)     # continuous problem with periodic BC
eig_disc_anal = np.zeros(N)     # discrete problem with periodic BC

for j in range(1, N+1):
    if j % 2 == 0:  # even eigenvalue
        arg = j
    else:           # odd eigenvalue
        arg = j-1
    eig_cont_anal[j-1] = arg**2
    eig_disc_anal[j-1] = np.sin(np.pi*arg/(2*N))**2

# scaling eigenvalues according to the corresponding analytical expressions
factor_cont = -(np.pi/L)**2
factor_disc = -(4./h**2)
eig_cont_anal *= factor_cont
eig_disc_anal *= factor_disc

# will calculate the relative error between the eigenvalues (except first lambda_0=0) obtained from
# a) analytic eigenvalues for the continuous  problem
# b) numerical eigenvalues for discrete problem with N points

error = np.zeros(N-1)
for j in range(1, N):
    error[j-1] = abs(A_eig[j]-eig_cont_anal[j])/abs(eig_cont_anal[j])

# print(eig_cont_anal); print(A_eig)
K = 4  # number of eigenvalues to print out a relative error
print(error[0:K])

其它你可能感兴趣的问题

上一篇了解设计模式下一篇python中BVP的哪些求解器是最好的？有没有比 scipy.integrate.solve_bvp 更好的东西？