利用扩散方程的矩阵求逆来计算沿杆的温度分布

Question

我的任务是随着时间的推移计算沿着1d杆的温度分布。我使用了this matrix，后面的时间等效公式是this。我的代码如下 - 我生成了一个节点网格，T似乎在每个位置都被计算：

import numpy as np
from scipy import linalg as lin
import matplotlib.pyplot as plt

#create 1d grid
N=50
xmax=50
x=np.linspace(0,xmax,N+1)
h=xmax/(N+1)
alpha=0.01
dt=1
print(x)
#create matrix

M=np.eye(N+1)
M=M*(1+(2*alpha*dt)/h**2)



for i in range(2,N+1):

     M[i-1,i]=-(alpha*dt)/h**2

for i in range(0,N-1):
     M[i+1,i]=-(alpha*dt)/h**2

M[0,0]=1
M[N-2,N-2]=1
print(M)


b=np.zeros((N+1,1))

    #set b.c.s
b[:,:]=20
b[N,:]=1000
print(b)

t=20




for j in range(0,t,1):

    LU, P = lin.lu_factor(M) 
    vector=lin.lu_solve((LU,P),b) 
    b=vector
    print(vector)


    plt.plot(x,vector, label='x')
#    plt.xlabel("k")
#    plt.ylabel("vector element")
#    pylab.legend(loc='upper right')
plt.show()

现在，图形应该看起来像this ...但是对于不同的t（时间）值，我得到的图形会非常急剧地增加（几乎呈指数级），并且根本不会表示沿着杆的温度分布！（注意：temp位于y轴上，位置位于x）上。我不确定我的代码是否有问题，或者我没有正确设置边界/容差。任何人都可以看到问题可能是什么？

谢谢

Answer 1

您的图形看起来是指数的，因为您将结束边界条件设置为1000，如果将初始边界设置为1000，那么它看起来更接近您建议的图形：

#set b.c.s
b[:,:]=20
b[0,:]=1000

另外我觉得你在这里弄错了：

M[0,0]=1
M[N-2,N-2]=1
print(M)

我认为它应该是：

M[0,0]=1
M[N,N]=1
print(M)

因为这将最后一个值设置为1而不是第二个到最后一个对角线元素。

再次改变M实际上是最后一个元素。我还认为问题的一部分是我们给它的时间量。你给出的示例图表是t = 1e9秒。因此，如果我们只使用非常大的解决方案的图表10，那么它看起来与该示例非常相似。这是我的代码：

t = 100000



soln = np.zeros((N+1, t))

for j in range(0, t, 1):

    LU, P = lin.lu_factor(M)
    vector = lin.lu_solve((LU, P), b)
    b = vector
    flat_list = []
    for sublist in vector:
        for item in sublist:
            flat_list.append(item)
    soln[:, j] = flat_list


for k in range(0, t, int(t/10)):
    plt.plot(x, soln[:, k], label='x')
#    plt.xlabel("k")
#    plt.ylabel("vector element")
#    pylab.legend(loc='upper right')


plt.show()