2024 XTU PDE 11题 对定解问题编程实现差分格式

11.\ 对定解问题
\[\begin{cases}u_{tt}-u_{xx}=(t^2-x^2)\sin(xt),\quad0<x<1,\quad0<t\leq1,\\[1ex]u(x,0)=0,\quad u_t(x,0)=x,\quad0\leq x\leq1,\\[1ex]u(0,t)=0,\quad u(1,t)=\sin(t),\quad0\leq t\leq1,\end{cases}\]
该问题真解为\[u(x, t) = \sin(xt).\]编程实现差分格式
\[\begin{aligned}&\delta_{t}^{2}u_{i}^{k}-\delta_{x}^{2}u_{i}^{k}=(t_k^2-x_i^2)\sin(x_{i}t_{k}),\quad1\leq i\leq m-1,\quad1\leq k\leq n-1,\\&u_{i}^{0}=0,\quad u_{i}^{1}=\tau x_{i},\quad1\leq i\leq m-1,\\&u_{0}^{k}=0,\quad u_{m}^{k}=\sin(t_{k}),\quad0\leq k\leq n.\end{aligned}\]
观察误差的收敛阶。

import numpy as np
from scipy import sparse as sym
from numpy import sin
import matplotlib.pyplot as plt
import sympy as sy

Nx = 80
Nt = 80

x = np.array([0, 1])
t = np.array([0, 1])
X = np.linspace(x[0], x[1], num=Nx + 1)

T = np.linspace(t[0], t[1], num=Nt + 1)
hx = X[1] - X[0]
tau = T[1] - T[0]
tmp = np.ones(Nx - 1)
tmp1 = np.ones(Nx - 2)
A1 = sym.diags([2 * (1 / tau ** 2 - 1 / hx ** 2) * tmp], [0], format='csc', dtype=np.float64)
A2 = sym.diags([1 / hx ** 2 * tmp1, 1 / hx ** 2 * tmp1], [-1, 1], format='csc', dtype=np.float64)
A = (A1 + A2) * tau ** 2

U = np.zeros((len(T), len(X)))

U[1, :] = tau * X
U[:, len(X) - 1] = sin(T)

F = np.empty_like(U)
x, t = sy.symbols('x t')
XX, TT = np.meshgrid(X, T)
U1 = np.zeros_like(U)


def f(x, t):
    return (t ** 2 - x ** 2) * sin(x * t)


F = f(XX, TT) * (tau ** 2)

for i in range(Nx - 1):
    partial_u = 1 / hx ** 2 * U[i + 1, -1]
    F[i + 1, -2] = F[i + 1, -2] + partial_u * tau ** 2


def solution(x, t):
    return sin(x * t)


U1 = solution(TT, XX)

for i in range(Nx - 1):
    U[i + 2, 1:-1] = A * U[i + 1, 1:-1] - U[i, 1:-1] + F[i + 1, 1:-1]

# print(U)
error = abs(U - U1)
e = np.max(np.max(error))
print(e)
plt.figure(figsize=(9, 9))
ax = plt.subplot(111, projection='3d')

ax.plot_surface(TT, XX, error,cmap='OrRd')
ax.set_xlabel("x")
ax.set_ylabel("t")
plt.show()