初学练习，看b站课程，教学为matlab代码，自己写的Python代码，后面会放b站课程链接，感兴趣的同学可以学习观看。

说明：Python初学者，代码可能不够漂亮，欢迎大家批评指正。本系列代码用notebook写的，有些图片显示的命令在pycharm中可能会报错，注释掉就行。

理论部分是课程笔记，有不明白的地方可以看b站课程（一位宝藏老师）。
数值方法3:偏微分方程1使用有限差分法解一维热传导（扩散）方程_哔哩哔哩_bilibilihttps://www.bilibili.com/video/BV1C7411f7fX?spm_id_from=333.999.0.0数值方法3:偏微分方程1使用有限差分法解一维热传导（扩散）方程(续）_哔哩哔哩_bilibili编程实现有限差分法解一维热传导（扩散）方程的例子。https://www.bilibili.com/video/BV1i7411f7Tr?spm_id_from=333.999.0.0数值方法3:偏微分方程1:使用有限差分法解一维热传导（扩散）方程3_哔哩哔哩_bilibili这次课讲存在热源（扩散源）的情况。并且讨论了不同的线性边界条件的处理。https://www.bilibili.com/video/BV1NE411w7a8?spm_id_from=333.999.0.0

 无热源

#!/usr/bin/env python
# coding: utf-8

# 数值方法3：偏微分方程1 使用有限差分法解一维热传导（扩散）方程
# 无热源情况

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
from matplotlib.animation import FFMpegWriter

# 参数
a = 1
dx = 0.02
X = 1
x = []
x = np.linspace(0,X,int(X/dx+1))
dt = 0.0001
T = 0.5
t = []
t = np.linspace(0,T,int(T/dt+1))

# 算子
A = (-2*np.eye(len(x),len(x),dtype=int)
     +np.diag(np.diag(np.ones([len(x),len(x)]),k=1),k=1)
     +np.diag(np.diag(np.ones([len(x),len(x)]),k=1),k=-1))

# 初始化
u = np.sin(np.pi*x)
# u = np.exp(-20*(np.power(x-0.5,2))) # 初始条件：高斯分布
# 边界
m1 = 0+0.1*np.sin(t)
m2 = 3-0.2*np.sin(8*t)

# 解方程组
u = np.mat(np.sin(np.pi*x)).T
# print(u.shape)
for n in range(len(t)-1):
    u = np.append(u,u[:,n]+a*a*dt/dx/dx*np.mat(A)*u[:,n],axis=1)
    u[0,n+1] = m1[n+1]
    u[-1,n+1] = m2[n+1]

# 绘制边界线
Y = np.array(u)
X = np.array(x)
plt.plot(X,Y[:,-1])
plt.title('FinalStates Temperature')
plt.savefig("FinalStatesT.png")  # 保存图片
plt.show()

# 3D曲面
get_ipython().run_line_magic('matplotlib', 'notebook')
# 准备数据
x_x, y_y = np.meshgrid(t,x)

# 绘制图片
fig = plt.figure("3D Surface", facecolor="lightgray")
plt.title("Temperature3D", fontsize=18)

# 设置为3D图片类型
ax3d = Axes3D(fig)

ax3d.set_xlabel("X")
ax3d.set_ylabel("Y")
ax3d.set_zlabel("Z")
plt.tick_params(labelsize=10)

ax3d.plot_surface(x_x, y_y, u, cstride=20, rstride=20, cmap="jet")

plt.savefig("Temperature3D.png")  # 保存图片
plt.show()

#初始化画布
fig = plt.figure()
# plt.xlim(0,X)
plt.ylim(0,np.max(u)) # np.max(psi)
plt.grid(ls='--')

ys = u[:,0]
Figure = plt.plot(x,ys,c='blue',alpha=0.8)[0]
#更新函数
def updata(num):
    ys = u[:,num]
    Figure.set_data(x,ys)
    return Figure

# animation库绘制动图
ani = animation.FuncAnimation(fig=fig,func=updata,frames=np.arange(0,int(len(t)/2)),interval=5) # frames帧数，interval间隔（帧）
plt.show()

## 保存为mp4，运行速度较慢，不保存时注释掉。
# mywriter = FFMpegWriter(fps=60)
# ani.save('TemperatureHistory.MP4',writer=mywriter)

 

有热源 

#!/usr/bin/env python
# coding: utf-8

# 数值方法3：偏微分方程1 使用有限差分法解一维热传导（扩散）方程
# 有热源情况

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
from matplotlib.animation import FFMpegWriter

# 参数
a = 1

dx = 0.02
X = 1
x = []
x = np.linspace(0,X,int(X/dx+1))

dt = 0.0001
T = 0.5
t = []
t = np.linspace(0,T,int(T/dt+1))

# 算子
A = (-2*np.eye(len(x),len(x),dtype=int)
     +np.diag(np.diag(np.ones([len(x),len(x)]),k=1),k=1)
     +np.diag(np.diag(np.ones([len(x),len(x)]),k=1),k=-1))

# 初始化
u = 0*x
# 边界
m1 = 0+0.0*np.sin(t)
m2 = 0-0.0*np.sin(8*t)
# 热源
f = np.mat(5*np.exp(-20*(np.power(x-0.5,2)))).T

# 解方程组
u = np.mat(np.sin(np.pi*x)).T
# print(u.shape)
for n in range(len(t)-1):
    u = np.append(u,np.add(u[:,n],np.add(a*a/dx/dx*np.mat(A)*u[:,n],f)*dt),axis=1)
    u[0,n+1] = m1[n+1]
    u[-1,n+1] = m2[n+1]

# 绘制边界线
Y = np.array(u)
X = np.array(x)
plt.plot(X,Y[:,-1])
plt.title('FinalStates Temperature')
plt.savefig("FinalStatesT.png")  # 保存图片
plt.show()

# 3D曲面
get_ipython().run_line_magic('matplotlib', 'notebook')
# 准备数据
x_x, y_y = np.meshgrid(t,x)

# 绘制图片
fig = plt.figure("3D Surface", facecolor="lightgray")
plt.title("Temperature3D", fontsize=18)

# 设置为3D图片类型
ax3d = Axes3D(fig)

ax3d.set_xlabel("X")
ax3d.set_ylabel("Y")
ax3d.set_zlabel("Z")
plt.tick_params(labelsize=10)

ax3d.plot_surface(x_x, y_y, u, cstride=20, rstride=20, cmap="jet")

plt.savefig("Temperature3D.png")  # 保存图片
plt.show()

#初始化画布
fig = plt.figure()
# plt.xlim(0,X)
plt.ylim(0,np.max(u)) # np.max(psi)
plt.grid(ls='--')

ys = u[:,0]
Figure = plt.plot(x,ys,c='blue',alpha=0.8)[0]
#更新函数
def updata(num):
    ys = u[:,num]
    Figure.set_data(x,ys)
    return Figure

# animation库绘制动图
ani = animation.FuncAnimation(fig=fig,func=updata,frames=np.arange(0,int(len(t)/2)),interval=2) # frames帧数，interval间隔（帧）
plt.show()

## 保存为mp4，运行速度较慢，不保存时注释掉。
# mywriter = FFMpegWriter(fps=60)
# ani.save('TemperatureHistory.MP4',writer=mywriter)

 

 

 第二类边界条件

#!/usr/bin/env python
# coding: utf-8

# 数值方法3：偏微分方程1 使用有限差分法解一维热传导（扩散）方程
# 第二类边界条件
#     两端绝热，图像在端点处与x轴平行，斜率为0
#         有热源情况，温度整体上升
#         无热源，温度逐渐均匀
#     一端绝热，一端有热流流入（斜率不为零）

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
from matplotlib.animation import FFMpegWriter

# 参数
a = 1

dx = 0.02
X = 1
x = []
x = np.linspace(0,X,int(X/dx+1))

dt = 0.0001
T = 0.1
t = []
t = np.linspace(0,T,int(T/dt+1))

# 算子
A = (-2*np.eye(len(x),len(x),dtype=int)
     +np.diag(np.diag(np.ones([len(x),len(x)]),k=1),k=1)
     +np.diag(np.diag(np.ones([len(x),len(x)]),k=1),k=-1))

# 初始化
u = 0*x
# 边界
# m1 = 0+0.0*np.sin(t)
m1 = -1+0.0*np.sin(t) # 热流流入
m2 = 0-0.0*np.sin(8*t)
# 热源
f = 10*np.mat(5*np.exp(-20*(np.power(x-0.5,2)))).T

# 解方程组
u = np.mat(np.sin(np.pi*x)).T
# print(u.shape)
for n in range(len(t)-1):
    u = np.append(u,np.add(u[:,n],np.add(a*a/dx/dx*np.mat(A)*u[:,n],f)*dt),axis=1)
    u[0,n+1] = u[1,n+1]-m1[n+1]*dx
    u[-1,n+1] = u[-2,n+1]+m2[n+1]*dx

get_ipython().run_line_magic('matplotlib', 'notebook # anaconda的jupyter notebook绘图界面，pycharm要注释掉')
# 绘制边界线
plt.figure(1) 
Y = np.array(u)
X = np.array(x)
plt.plot(X,Y[:,-1])
plt.title('FinalStates Temperature')
plt.savefig("FinalStatesT.png")  # 保存图片
plt.show()

# 3D曲面
plt.figure(2) 
# 准备数据
x_x, y_y = np.meshgrid(t,x)

# 绘制图片
fig = plt.figure("3D Surface", facecolor="lightgray")
plt.title("Temperature3D", fontsize=18)

# 设置为3D图片类型
ax3d = Axes3D(fig)

ax3d.set_xlabel("X")
ax3d.set_ylabel("Y")
ax3d.set_zlabel("Z")
plt.tick_params(labelsize=10)

ax3d.plot_surface(x_x, y_y, u, cstride=20, rstride=20, cmap="jet")

plt.savefig("Temperature3D.png")  # 保存图片
plt.show()

#初始化画布
fig = plt.figure(3)
# plt.xlim(0,X)
plt.ylim(0,np.max(u)) # np.max(psi)
plt.grid(ls='--')

ys = u[:,0]
Figure = plt.plot(x,ys,c='blue',alpha=0.8)[0]
#更新函数
def updata(num):
    ys = u[:,num]
    Figure.set_data(x,ys)
    return Figure

# animation库绘制动图
ani = animation.FuncAnimation(fig=fig,func=updata,frames=np.arange(0,int(len(t)/2)),interval=5) # frames帧数，interval间隔（帧）
plt.show()

## 保存为mp4，运行速度较慢，不保存时注释掉
# mywriter = FFMpegWriter(fps=20)
# ani.save('TemperatureHistory.MP4',writer=mywriter)

 

 第三类边界条件

#!/usr/bin/env python
# coding: utf-8

# 数值方法3：偏微分方程1 使用有限差分法解一维热传导（扩散）方程
# 第三类边界条件

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
from matplotlib.animation import FFMpegWriter

# 参数
a = 1
dx = 0.02
X = 1
x = []
x = np.linspace(0,X,int(X/dx+1))
dt = 0.0001
T = 0.1
t = []
t = np.linspace(0,T,int(T/dt+1))
c = 1

# 算子
A = (-2*np.eye(len(x),len(x),dtype=int)
     +np.diag(np.diag(np.ones([len(x),len(x)]),k=1),k=1)
     +np.diag(np.diag(np.ones([len(x),len(x)]),k=1),k=-1))

# 初始化
u = 0*x
# 边界
m1 = 0+0.0*np.sin(t)
# m1 = -1+0.0*np.sin(t) # 热流流入
m2 = 0-0.0*np.sin(8*t)
# 热源
f = 10*np.mat(5*np.exp(-20*(np.power(x-0.5,2)))).T

# 解方程组
u = np.mat(np.sin(np.pi*x)).T
# print(u.shape)
for n in range(len(t)-1):
    u = np.append(u,np.add(u[:,n],np.add(a*a/dx/dx*np.mat(A)*u[:,n],f)*dt),axis=1)
    u[0,n+1] = (c*u[1,n+1]-m1[n+1]*dx)/(c-dx)
    u[-1,n+1] = (1-dx/c)*u[-2,n+1]+m2[n+1]*dx/c


# anaconda的jupyter notebook绘图界面，pycharm要注释掉
get_ipython().run_line_magic('matplotlib', 'notebook')

# 绘制边界线
plt.figure(1) 
Y = np.array(u)
X = np.array(x)
plt.plot(X,Y[:,-1])
plt.title('FinalStates Temperature')
plt.savefig("FinalStatesT.png")  # 保存图片
plt.show()

# 3D曲面
# 准备数据
x_x, y_y = np.meshgrid(t,x)

# 绘制图片
fig2 = plt.figure("3D Surface", facecolor="lightgray")
plt.title("Temperature3D", fontsize=18)

# 设置为3D图片类型
ax3d = Axes3D(fig2)

ax3d.set_xlabel("X")
ax3d.set_ylabel("Y")
ax3d.set_zlabel("Z")
plt.tick_params(labelsize=10)

ax3d.plot_surface(x_x, y_y, u, cstride=20, rstride=20, cmap="jet")

plt.savefig("Temperature3D.png")  # 保存图片
plt.show()

#初始化画布
fig3 = plt.figure(3)
# plt.xlim(0,X)
plt.ylim(0,np.max(u)) # np.max(psi)
plt.grid(ls='--')

ys = u[:,0]
Figure = plt.plot(x,ys,c='blue',alpha=0.8)[0]
#更新函数
def updata(num):
    ys = u[:,num]
    Figure.set_data(x,ys)
    return Figure

# animation库绘制动图
ani = animation.FuncAnimation(fig=fig3,func=updata,frames=np.arange(0,int(len(t)/2)),interval=5) # frames帧数，interval间隔（帧），越小速度越快
plt.show()

## 保存为mp4，运行速度较慢，不保存时注释掉。fps越小速度越慢
# mywriter = FFMpegWriter(fps=20)
# ani.save('TemperatureHistory.MP4',writer=mywriter)

 

# python
Logo
华为云开发者联盟

为开发者提供学习成长、分享交流、生态实践、资源工具等服务，帮助开发者快速成长。

更多推荐

cover

技术解密Java Chassis 3超实用的可观测性

avatar 华为云开发者联盟
cover

详解Spring Boot的RedisAutoConfiguration配置

avatar 华为云开发者联盟
cover

揭秘华为云运维中心，如何守护全球10亿用户的智慧生活体验？

avatar 华为云开发者联盟