# -*- coding: utf-8 -*-
"""
Created on Sat Jul 6 08:18:27 2024
@author: chlj
"""
# Solving partial differential equations
# 偏微分方程数值解法
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
# 抛物线方程
# U_t = D * U_xx
# 隐式方法边界条件U(0,t) = 0, U(1,t) = 0
# 初始条件函数 U(x,0)
def funcUx0(x):
u0 = (np.sin(2 * np.pi * x)) ** 2
return u0
M = 101 #空间步数
N = 101 #时间步数
h = 1/(M-1) #x的步长
k = 1/(N-1) #t的步长
D = 1 #扩散系数
sigma = D*k /(h*h) #σ
#设置矩阵
A = np.zeros((M-2, M-2))
w = np.zeros((M-2, N))
for i in range(0,M-2):
for j in range(0,M-2):
if j == i:
A[i][j] = 1 + 2 * sigma
elif j == i+1:
A[i][j] = -sigma
elif j == i-1:
A[i][j] = -sigma
'''
[[ 201. -100. 0. ... 0. 0. 0.]
[-100. 201. -100. ... 0. 0. 0.]
[ 0. -100. 201. ... 0. 0. 0.]
...
[ 0. 0. 0. ... 201. -100. 0.]
[ 0. 0. 0. ... -100. 201. -100.]
[ 0. 0. 0. ... 0. -100. 201.]]
'''
A_inv = np.linalg.inv(A) #求逆
b = np.zeros((M-2, 1))
for i in range(1,M-1):
b[i-1][0] = funcUx0(k*i)
w[:,0]= b.T
for i in range(0,N-2):
b = np.dot(A_inv,b)
w[:,i+1] = b.T
c = np.zeros((1,N)) #边界条件
w = np.r_[c,w,c] #添加边界条件
#画网格
x, t = np.meshgrid(np.linspace(0, 1, M),
np.linspace(0, 1, N))
#画图
f = plt.figure('Wireframe', facecolor='lightgray')
ax1 = f.add_subplot(111,projection='3d')
ax1.set_xlabel('X', fontsize=14)
ax1.set_ylabel('T', fontsize=14)
ax1.set_zlabel('Z', fontsize=14)
ax1.set_zlim(-1, 1)
ax1.plot_wireframe(x, t, w, linewidth=0.5)
ax1.view_init(elev=20, azim=-120)
plt.savefig('抛物线方程.jpg',dpi = 200)
plt.show()
- 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.
- 38.
- 39.
- 40.
- 41.
- 42.
- 43.
- 44.
- 45.
- 46.
- 47.
- 48.
- 49.
- 50.
- 51.
- 52.
- 53.
- 54.
- 55.
- 56.
- 57.
- 58.
- 59.
- 60.
- 61.
- 62.
- 63.
- 64.
- 65.
- 66.
- 67.
- 68.
- 69.
- 70.
- 71.
- 72.
- 73.
- 74.
- 75.
- 76.
- 77.
- 78.
- 79.
- 80.
- 81.
- 82.