Python——求解一维热传导

1、问题

     求一维热传导方程混合问题：

     求解其的数值解，取N=10，h=0.1，计算到K=36为止.

 2、程序

#python
import numpy as np
import matplotlib.pyplot as plt
def f(a,x1,x2,x3):
    y=a*x3+(1-2*a)*x2+a*x1
    return y
h=0.1;a=1/6;p=1/600
x=np.zeros((37,6));
for i in range(6):
    x[0,i]=4*i*h*(1-i*h)
for j in range(1,37):
    for n in range(1,6):
        if n+1!=6:
            x[j,n]=f(a,x[j-1,n-1],x[j-1,n],x[j-1,n+1])
        else:
            x[j,n]=f(a,x[j-1,n-1],x[j-1,n],x[j-1,n-1])
y1=[0,1,2,3,4,5]
t=np.array([])
for o in range(37):
    t1=[p*o,p*o,p*o,p*o,p*o,p*o]
    t=np.r_[t,t1]
t=np.array(t.reshape(37,-1))
po=['position:0','position:1','position:2','position:3','position:4','position:5']
for by in range(6):
    x2=x[:,by]
    y2=t[:,by]
    plt.plot(y2,x2,label=po[by])
plt.xlabel('Time')
plt.ylabel('Heat')
plt.title('Heat changes over time')
plt.legend()
plt.show()

3、结果

运行结果如Fig1.

Fig1. Position(0~5)表示空间位置，y-axis 表示热量/焦耳，x-axis表示时间/秒

 

4、实验总结

（1）本次实验解决偏微分方程所用到的方法是差分法，核心思想就是将一阶差商和二阶差商代替一阶偏分和二阶偏分；

（2）本次实验程序的实现并不困难，掌握核心方法，利用边值条件，对不同空间位置的点进行迭代就可以较为容易的出现结果；

（3）从设计的程序运行结果来看，随着时间的增大，热量在不断的减少，这种规律是符号我们经验的，具有一定的可信性。