python有限差分法求解一维热传导方程

​1、方程及其离散

1.1一维热传导方程

\left\{\begin{matrix} \frac{\partial u\left(x,t \right )}{\partial x}=\frac{\partial u\left(x,t \right )}{\partial x^2}\\ u\left(x,0 \right )=4x\left(3-x \right ) \\ u\left(0,t \right )=0 \\ u\left(3,t \right )=0 \end{matrix}\right.

1\leqslant t\leqslant 1000,0\leqslant x\leqslant 1

1.2离散化

设定步长h=0.1,\tau =0.0001,依据上述方程得到递推关系:

2、python求解实现

import numpy as np
import matplotlib.pyplot as plt

h = 0.1#空间步长
N =30#空间步数
dt = 0.0001#时间步长
M = 10000#时间的步数
A = dt/(h**2) #lambda*tau/h^2
U = np.zeros([N+1,M+1])#建立二维空数组
Space = np.arange(0,(N+1)*h,h)#建立空间等差数列,从0到3,公差是h

#边界条件
for k in np.arange(0,M+1):
    U[0,k] = 0.0
    U[N,k] = 0.0

#初始条件
for i in np.arange(0,N):
    U[i,0]=4*i*h*(3-i*h)

#递推关系
for k in np.arange(0,M):
    for i in np.arange(1,N):
        U[i,k+1]=A*U[i+1,k]+(1-2*A)*U[i,k]+A*U[i-1,k]
#不同时刻的温度随空间坐标的变化
plt.plot(Space,U[:,0], 'g-', label='t=0',linewidth=1.0)
plt.plot(Space,U[:,3000], 'b-', label='t=3/10',linewidth=1.0)
plt.plot(Space,U[:,6000], 'k-', label='t=6/10',linewidth=1.0)
plt.plot(Space,U[:,9000], 'r-', label='t=9/10',linewidth=1.0)
plt.plot(Space,U[:,10000], 'y-', label='t=1',linewidth=1.0)
plt.ylabel('u(x,t)', fontsize=20)
plt.xlabel('x', fontsize=20)
plt.xlim(0,3)
plt.ylim(-2,10)
plt.legend(loc='upper right')
plt.show()

  • 0
    点赞
  • 16
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
一维热传导方程可以表示为: $$\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$$ 其中,$u(x,t)$ 表示温度分布,$\alpha$ 为热扩散系数。 使用有限差分方法求解一维热传导方程可以分为以下几个步骤: 1. 离散化:将时间和空间分别离散化为 $t_i = i \Delta t$ 和 $x_j = j \Delta x$,其中 $\Delta t$ 和 $\Delta x$ 分别表示时间和空间的步长。 2. 近似:将偏导数用差分来近似,例如: $$\frac{\partial u}{\partial t} \approx \frac{u_{j,i+1} - u_{j,i}}{\Delta t}$$ $$\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{j+1,i} - 2u_{j,i} + u_{j-1,i}}{\Delta x^2}$$ 3. 替换:将近似后的偏导数代入原方程,得到差分方程: $$u_{j,i+1} = u_{j,i} + \alpha \frac{\Delta t}{\Delta x^2} (u_{j+1,i} - 2u_{j,i} + u_{j-1,i})$$ 4. 边界条件:确定边界条件,例如固定边界条件 $u(0,t) = u(L,t) = 0$。 5. 时间推进:使用差分方程逐步推进时间,例如使用显式欧拉法: $$u_{j,i+1} = u_{j,i} + \alpha \frac{\Delta t}{\Delta x^2} (u_{j+1,i} - 2u_{j,i} + u_{j-1,i})$$ 6. 初始条件:确定初始条件,例如 $u(x,0) = f(x)$。 下面是一个 Python 实现的示例代码: ```python import numpy as np import matplotlib.pyplot as plt # 设置参数 L = 1.0 # 空间长度 T = 0.1 # 时间长度 alpha = 1.0 # 热扩散系数 N = 100 # 空间网格数 M = 100 # 时间网格数 dx = L / N # 空间步长 dt = T / M # 时间步长 # 初始化温度分布 u = np.zeros((N+1, M+1)) u[0,:] = u[N,:] = 0.0 u[:,0] = np.sin(np.pi*np.arange(N+1)/(N+1)) # u[:,0] = np.exp(-100*(np.arange(N+1)/(N+1) - 0.5)**2) # 高斯分布 # 使用有限差分方法求解热传导方程 for i in range(M): for j in range(1, N): u[j,i+1] = u[j,i] + alpha*dt/dx**2 * (u[j+1,i] - 2*u[j,i] + u[j-1,i]) # 绘制温度分布曲线 x = np.linspace(0, L, N+1) t = np.linspace(0, T, M+1) X, T = np.meshgrid(x, t) fig = plt.figure(figsize=(8,6)) ax = fig.add_subplot(1,1,1,projection='3d') ax.plot_surface(X, T, u.T, cmap='coolwarm') ax.set_xlabel('x') ax.set_ylabel('t') ax.set_zlabel('u') plt.show() ``` 其中,初始化温度分布的代码可以根据实际情况进行修改。绘制温度分布曲线的代码使用了 `matplotlib` 库中的 `plot_surface` 函数。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值