1
%reset -f

解偏微分方程——差分法

$$\dfrac{\partial u}{\partial t}=\nu\dfrac{\partial^2u}{\partial x^2}+\nu\dfrac{\partial^2u}{\partial y^2}$$
$$\dfrac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t}=\nu\dfrac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}+\nu\dfrac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}$$
$$u_{i,j}^{n+1}=u_{i,j}^n+\dfrac{\nu\Delta t}{\Delta x^2}\big(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n\big)+\dfrac{\nu\Delta t}{\Delta y^2}\big(u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n\big)$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

nx = 100  # 离散化x方向格数
ny = 100  # y方向格数
nu = 1.  # 
dx = 2 / (nx - 1)  # 离散的deltax 即网格长
dy = 2 / (ny - 1)  # 离散的deltay 网格宽
t = 0.1  # 

x = np.linspace(-1,1,nx)
y = np.linspace(-1,1,ny)


x1 = x.reshape((-1, 1))
y1 = y.reshape((1, -1))
tmpx1 = np.sin(np.pi*x1)
tmpy1 = np.sin(np.pi*y1)
u = np.matmul(tmpx1, tmpy1)


u[0, :] = 0
u[-1, :] = 0
u[:, 0] = 0
u[:, -1] = 0


fig = plt.figure()
ax = fig.gca(projection='3d')
xa, ya = np.meshgrid(x, y)
surf2 = ax.plot_surface(xa,ya,u[:], cmap='viridis')
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
C:\Users\w1858\AppData\Local\Temp\ipykernel_11756\328840734.py:30: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax = fig.gca(projection='3d')





Text(0.5, 0.5, '$y$')

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
nt= 1000
dt = t / nt
# 定义函数进行求解
for n in range(nt + 1):
	un = u.copy()
	u[1:-1,1:-1] = (un[1:-1,1:-1] +
		             ((nu * dt) / (dx**2)) * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]) +
		             ((nu * dt) / (dy**2)) * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) )
	u[0, :] = 0
	u[-1, :] = 0
	u[:, 0] = 0
	u[:, -1] = 0

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xa,ya,u[:], rstride=1, cstride=1, cmap='viridis',
        linewidth=0, antialiased=False)
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-0.2, 0.2)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.show()
C:\Users\w1858\AppData\Local\Temp\ipykernel_11756\1993782575.py:15: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax = fig.gca(projection='3d')

png