2024 湘潭大学 PDE 第四题五点差分格式_五点差分求解pde-CSDN博客

利用五点差分格式求解 $\left\{ \begin{aligned} -\Delta u+u=& f(x,y),\quad(x,y)\in\Omega=(0,2)\times(0,1), \\ u|_{\partial\Omega}=& g(x,y). \end{aligned} \right.$

该问题的真解为 $u(x,y)=e^{x}\sin(2\pi y),f(x,y)$ 由真解 $u(x,y)$ 代入方程计算得到。分别取步长 $(h_{x},h_{y})=(\frac{1}{10},\frac{1}{20}),(\frac{1}{20},\frac{1}{40}),(\frac{1}{40},\frac{1}{80}),(\frac{1}{80},\frac{1}{160})$ 进行计算，观察误差 $\|e\|_C\text{、}\|e\|_0$ 的收敛阶，画出真解曲面图、数值解曲面图和误差曲面图。

import numpy as np
from scipy import sparse as sym
from numpy import sin, cos, pi, exp
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt


class CDModel:
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def space_grid(self, Nx, Ny):
        X = np.linspace(self.x[0], self.x[1], Nx + 1)
        Y = np.linspace(self.y[0], self.y[1], Ny + 1)
        hx = (self.x[1] - self.x[0]) / Nx
        hy = (self.y[1] - self.y[0]) / Ny
        return X, Y, hx, hy

    def f(self, X, Y):
        return 4 * pi ** 2 * exp(X) * sin(2 * pi * Y)

    def solution(self, X, Y):
        return exp(X) * sin(2 * pi * Y)

    def partial_u(self, Nx, Ny):
        """有点错误，这个需要自己去看边界，名称不对"""
        X, Y, hx, hy = self.space_grid(Nx, Ny)
        ub = sin(2 * pi * Y)
        ut = exp(2) * sin(2 * pi * Y)
        ul = np.zeros(len(X))
        ur = np.zeros(len(X))
        return ul, ur, ut, ub

    def source(self, Nx, Nt):
        X, Y, hx, hy = self.space_grid(Nx, Nt)
        XX, YY = np.meshgrid(X, Y)
        f_1 = self.f(XX, YY)
        return f_1

    def solution_matrix(self, Nx, Nt):
        X1, Y, hx, hy = self.space_grid(Nx, Nt)
        XX, YY = np.meshgrid(X1, Y)
        U = self.solution(XX, YY)
        return U

    def cd_Iterative(self, Nx, Ny):
        X, Y, hx, hy = self.space_grid(Nx, Ny)
        tmp = np.ones(Nx - 1)
        tmp1 = np.ones(Nx - 2)
        A = sym.diags([(1 + 2 / hx ** 2 + 2 / hy ** 2) * tmp, (-1 / hx ** 2) * tmp, (-1 / hx ** 2) * tmp], [0, 1, -1],
                      format='csc', dtype=np.float64)
        E = sym.eye(Nx - 1, format='csc')
        B = sym.diags(-1 / hy ** 2 * tmp, 0, format='csc', dtype=np.float64)
        E_1 = sym.diags([1 * tmp1, 1 * tmp1], [-1, 1], format='csc', dtype=np.float64)
        A_1 = sym.kron(E_1, B) + sym.kron(E, A)
        ul, ur, ut, ub = self.partial_u(Nx, Ny)
        f = self.source(Nx, Ny)
        f[:, 1] += 1 / hx ** 2 * ub
        f[:, -2] += 1 / hx ** 2 * ut
        f = f[1:-1, 1:-1].flatten()
        u = np.zeros((len(X), len(Y)))
        u[1:-1, 1:-1] = spsolve(A_1, f).reshape((len(X) - 2, len(Y) - 2))
        u[-1, :] = ul
        u[0, :] = ur
        u[:, 0] = ub
        u[:, -1] = ut
        return u

    def error(self, Nx, Nt):
        u_1 = self.solution_matrix(Nx, Nt)
        u = self.cd_Iterative(Nx, Nt)
        error = u_1 - u
        return error


def text():
    x = np.array([0, 2])
    y = np.array([0, 1])
    data = CDModel(x, y)
    u = data.cd_Iterative(Nx, Ny)
    Error = data.error(Nx, Ny)
    Error = abs(Error)
    U = data.solution_matrix(Nx, Ny)
    return u, Error, U


Nx = 320
Ny = 320
hx = 2 / Nx
hy = 1 / Ny

u1, Error, U = text()
ofanshu = np.sqrt(sum(np.multiply(Error ** 2, hx * hy)))

a = np.max(np.max(Error))
print(a)
X = np.linspace(0, 2, Nx + 1)
Y = np.linspace(0, 1, Ny + 1)
plt.figure(figsize=(15, 9))
ax = plt.subplot(121, projection='3d')
XX, YY = np.meshgrid(X, Y)
ax.plot_surface(XX, YY, U, cmap='hot')
ax1 = plt.subplot(122, projection='3d')
ax1.plot_surface(XX,YY,Error)
plt.show()