2024 湘潭大学 PDE数值解 中心差分格式 第三题

应用中心差分格式计算如下定解问题：

该问题的真解为由真解$u(x)$代入方程计算得到。

分别取步长进行计算，观察误差的收敛阶，并画出真解曲线、数值解曲线和误差曲线。

import numpy as np
import sympy as sy
import matplotlib.pyplot as plt
import matplotlib

x = sy.symbols('x')
u = sy.exp(-x ** 2) * (1 - x ** 2)
f = -sy.diff(u, x, 2) + u
h = 1 / 10
N = int(2 / h)
b = np.arange(-1, 1, h)
b1 = []
for i in range(N - 1):
    c = f.evalf(subs={x: b[i + 1]})
    b1.append(c)

u2=[]
for i in range(N - 1):
    c = u.evalf(subs={x: b[i + 1]})
    u2.append(c)

def tridiagonal_solver(a, b, c, d):
    n = len(a)
    alpha = [0] * n
    beta = [0] * n
    x = [0] * n

    # 初始化alpha和beta
    alpha[0] = b[0] / a[0]
    beta[0] = d[0] / a[0]

    # 迭代计算alpha和beta
    for i in range(1, n):
        alpha[i] = b[i] / (a[i] - c[i - 1] * alpha[i - 1])
        beta[i] = (d[i] - c[i - 1] * beta[i - 1]) / (a[i] - c[i - 1] * alpha[i - 1])

    # 回代求解x
    x[n - 1] = beta[n - 1]
    for i in range(n - 2, -1, -1):
        x[i] = beta[i] - alpha[i] * x[i + 1]

    return x


# 示例用例
a11 = np.array([1 + 2 / h ** 2] * (N - 1))
b11 = np.array([-1 / h ** 2] * (N - 1))
c11 = np.array([-1 / h ** 2] * (N - 1))
b1 = np.array(b1)
u1 = tridiagonal_solver(a11, b11, c11, b1)
u1 = np.array(u1)
e1=[]
for i in range(N-1):
    a = abs(u2[i]-u1[i])
    e1.append(a)

e1 = np.array(e1)


font = {'family': 'Microsoft YaHei', 'weight': 'bold', 'size': 12}
plt.rcParams['font.sans-serif'] = 'Microsoft YaHei'
plt.rcParams['axes.unicode_minus'] = False
matplotlib.rc('font', **font)
a = np.linspace(-1, 1, N-1)
x = np.array(a)
y = np.array(u1)


fig = plt.figure(figsize=(10, 4))
ax = fig.add_subplot(1, 1, 1)

h1, = plt.plot(x, y, color='b', label='解析解')
h2, = plt.plot(x, u2, color='r', label='真解')
ax.set_xlabel('x')
ax.set_ylabel('真解和解析解')

ax2 = ax.twinx()
h3 = plt.plot(x, e1, color='g', label='误差')
ax2.set_ylabel('误差')

ax2.spines['left'].set_color('r')
ax2.spines['right'].set_color('g')
ax.tick_params(axis='y', colors='r')
ax2.tick_params(axis='y', colors='g')
ax.yaxis.label.set_color('r')
ax2.yaxis.label.set_color('g')

ax.legend([h1, h2], ['解析解', '真解', '误差'])
plt.show()