利用五点差分格式求解
该问题的真解为由真解
代入方程计算得到。分别取步长
进行计算,观察误差
的收敛阶,画出真解曲面图、数值解曲面图和误差曲面图。
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()