Использование матричной инверсии уравнения диффузии для расчета распределения температуры вдоль стержня - PullRequest
Новая
       122

Использование матричной инверсии уравнения диффузии для расчета распределения температуры вдоль стержня

0 голосов
/ 04 марта 2019

Моя задача - вычислить распределение температуры вдоль стержня 1d с течением времени. Я использовал эту матрицу , и формула, эквивалентная времени в обратном направлении, равна это .Мой код выглядит следующим образом: я создал сетку узлов, и кажется, что 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, а позиция нах).Я не уверен, что с моим кодом что-то не так или я неправильно устанавливаю границы / допуски.Кто-нибудь может увидеть, в чем может быть проблема?

Спасибо

1 Ответ

0 голосов
/ 04 марта 2019

Ваш график выглядит экспоненциальным, потому что вы устанавливаете конечное граничное условие равным 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, а не вторым по последнему диагональному элементу.

После изменения М снова, чтобы фактически быть последним элементом.Я также думаю, что частью проблемы является количество времени, которое мы уделяем этому.Пример графика, который вы дали, уходит в t = 1e9 секунд.Таким образом, если мы сделаем t очень большим, чем только график 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()
...