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

初学练习，看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 = 1dx = 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 = 1dx = 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]*dxu[-1,n+1] = u[-2,n+1]+m2[n+1]*dxget_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)

数值方法3：偏微分方程1：使用有限差分法解一维热传导（扩散）方程相关推荐

一维热传导方程 matlab隐式解,一维热传导偏微分方程的数值解的matlab程序问题出在哪儿？...
我现在编写了一个求解一维热传导的偏微分方程,调程序都调了好多天了不知道问题在哪儿,求各位高手帮忙看一下好么? 我在此表示万分感谢需求解的方程看图片,我的程序如下 %---------------- ...
MATLAB模拟导热过程,一维热传导MATLAB模拟.doc
PAGE 昆明学院 2015 届毕业设计(论文) 设计(论文)题目一维热传导问题的数值解法及其MATLAB模拟子课题题目无姓名伍有超学号 201117030225 所属系物 ...
matlab解薛定谔方程,有限差分法解薛定谔方程与MATLAB实现
第30卷第3期高师理科学刊 Vol. 30 No.3 2010年 5月 Journal of Science of Teachers′ College and University Ma ...
有限差分法的一维扩散MATLAB,一维扩散方程的有限差分法matlab
用matlab编程实现一维扩散方程的有限差分法 1 一维扩散方程的有限差分法 --计算物理实验作业七陈万物理学2013级 130******** ● 题目: 编程求解一维扩散方程的解 ⎪⎪⎪⎪⎩ ...
利用有限元法求解一维热传导问题
我们都知道热传导方程是大名鼎鼎的傅里叶提出的,而对于此方程的解也可以用傅里叶变换的方法得到,但笔者已经记不起该怎么求解了(尽管上学期才刚学过数理方法啊).但是无伤大雅,因为我已经掌握了(其实还没有)对 ...
SIMPLE算法求解多孔介质的一维流动控制方程
SIMPLE算法求解多孔介质的一维流动控制方程问题介绍求解思想压力修正方法的基本思想两个关键问题求解步骤及说明疑惑: a e a_e ae表示动量方程离散系数?怎么求解? 本案例求解分析 ...
有限差分法求解不可压NS方程
网上关于有限差分法解NS方程的程序实现不尽完备,这里是一些补充注解现有的优秀资料理论向 [1]如何从物理意义上理解NS方程? - 知乎 [2]NS方程数值解法:投影法的简单应用 - 知乎 [3][ ...
热传递 matlab,一维热传导MATLAB模拟.pdf
昆明学院 2015 届毕业设计(论文) 设计(论文)题目一维热传导问题的数值解法及其 MATLAB 模拟子课题题目无姓名伍有超学号 201117030225 所属系物理科学 ...
用Matlab解二阶非齐次微分方程
用Matlab解二阶非齐次微分方程大纲函数代码大纲用Matlab解二阶非齐次微分方程,网上很多麻烦又累赘又无用的东西,一句话解决的事. 函数 dsolve('a','b','c'):解微分方 ...

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

数值方法3：偏微分方程1：使用有限差分法解一维热传导（扩散）方程相关推荐

最新文章

热门文章