再将精确解画出三维曲面图
时间: 2023-08-04 15:08:50 浏览: 50
要绘制精确解的三维曲面图,我们需要首先找到精确解的解析表达式,然后使用该表达式计算在每个离散点上的精确解值。在这个问题中,已经给出了精确解的表达式为 u(x, y) = e^xsin(πy)。
下面是修改后的代码,将精确解的三维曲面图与数值解进行对比:
```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def solve_by_elimination(m1, m2):
hx = 2 / m1
hy = 1 / m2
N = (m1 + 1) * (m2 + 1)
A = np.zeros((N, N))
b = np.zeros(N)
# 构建线性方程组系数矩阵和右侧向量
for i in range(1, m1):
for j in range(1, m2):
k = i * (m2 + 1) + j
A[k, k] = -2 / hx**2 - 2 / hy**2
A[k, k-1] = 1 / hx**2
A[k, k+1] = 1 / hx**2
A[k, k-(m2+1)] = 1 / hy**2
A[k, k+(m2+1)] = 1 / hy**2
b[k] = (np.pi**2 - 1) * np.exp(i * hx) * np.sin(np.pi * j * hy)
# 处理边界条件
for j in range(m2 + 1):
k = j
A[k, k] = 1
b[k] = np.sin(np.pi * j * hy)
k = m1 * (m2 + 1) + j
A[k, k] = 1
b[k] = np.exp(2) * np.sin(np.pi * j * hy)
for i in range(m1 + 1):
k = i * (m2 + 1)
A[k, k] = 1
b[k] = 0
k = i * (m2 + 1) + m2
A[k, k] = 1
b[k] = 0
# 求解线性方程组
u = np.linalg.solve(A, b)
return u.reshape((m1+1, m2+1))
def solve_by_gauss_seidel(m1, m2):
hx = 2 / m1
hy = 1 / m2
N = (m1 + 1) * (m2 + 1)
u = np.zeros(N)
# 构建线性方程组系数矩阵和右侧向量
A = np.zeros((N, N))
b = np.zeros(N)
for i in range(1, m1):
for j in range(1, m2):
k = i * (m2 + 1) + j
A[k, k] = -2 / hx**2 - 2 / hy**2
A[k, k-1] = 1 / hx**2
A[k, k+1] = 1 / hx**2
A[k, k-(m2+1)] = 1 / hy**2
A[k, k+(m2+1)] = 1 / hy**2
b[k] = (np.pi**2 - 1) * np.exp(i * hx) * np.sin(np.pi * j * hy)
# 迭代求解线性方程组
max_iter = 1000 # 最大迭代次数
tol = 0.5e-10 # 收敛精度
for iter in range(max_iter):
u_new = u.copy()
for i in range(1, m1):
for j in range(1, m2):
k = i * (m2 + 1) + j
u_new[k] = (b[k] - (A[k, :k] @ u_new[:k] + A[k, k+1:] @ u[k+1:])) / A[k, k]
if np.linalg.norm(u_new - u, ord=np.inf) < tol:
break
u = u_new
return u.reshape((m1+1, m2+1))
# 设置网格参数
m1 = 10 # x方向网格点数
m2 = 10 # y方向网格点数
# 使用消元法求解差分方程组
u_elimination = solve_by_elimination(m1, m2)
# 使用Gauss-Seidel迭代法求解差分方程组
u_gauss_seidel = solve_by_gauss_seidel(m1, m2)
# 计算精确解
x = np.linspace(0, 2, m1+1)
y = np.linspace(0, 1, m2+1)
X, Y = np.meshgrid(x, y)
u_exact = np.exp(X) * np.sin(np.pi * Y)
# 绘制三维曲面图
fig = plt.figure(figsize=(12, 4))
# 绘制消元法求解结果的三维曲面图
ax = fig.add_subplot(131, projection='3d')
ax.plot_surface(X, Y, u_elimination, cmap='hot')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u(x, y)')
ax.set_title('Numerical Solution by Elimination')
# 绘制Gauss-Seidel迭代法求解结果的三维曲面图
ax = fig.add_subplot(132, projection='3d')
ax.plot_surface(X, Y, u_gauss_seidel, cmap='hot')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u(x, y)')
ax.set_title('Numerical Solution by Gauss-Seidel')
# 绘制精确解的三维曲面图
ax = fig.add_subplot(133, projection='3d')
ax.plot_surface(X, Y, u_exact, cmap='hot')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u(x, y)')
ax.set_title('Exact Solution')
plt.tight_layout()
plt.show()
```
通过以上代码,我们可以同时绘制消元法求解结果、Gauss-Seidel迭代法求解结果和精确解的三维曲面图。这样可以直观地比较数值解和精确解之间的差异。希望这次的代码能满足你的需求,如果还有其他问题,请随时提问!