使用 python 求解热方程

Question

我正在尝试求解热方程，然后绘制如下图表：

import numpy as np
import matplotlib.pyplot as plt

L = 0.04     # soil depth
Te = 28    # ext tempe
Tn = 36.6      # tempe at x = L
D = 0.75    #  m2.s-1
tau = 232920  # total lengh of experiment
dt = 600   
dx = 0.001  
Tm=38   #average temperature
T=24*3600 #period
w=(2*pi)/T

Nt = int(tau/dt)+1  
Nx = int(L/dx)+1    

x = [j*dt for j in range(Nx)]       
T = np.zeros((Nt,Nx))               
for j in range(Nx):
    T[0,j] = Te    # T(x,t=0) = Te
T[0,0] = Te         # T(x=0,t=0) = Te
for i in range(1,Nt): 
    T[i,0] = Te+Tm*np.cos(w*i)     # T(x=0,t) = ...
    T[i,-1] = Tn    # T(x=L,t) = Tn 


for i in range(Nt-1):   
    for j in range(1,Nx-1):    
        T[i+1,j]=T[i,j]+D*dt/(dx*dx)*(T[i,j+1]+T[i,j-1]-2*T[i,j])

t=np.arange(0, 232920, 10*60)
plt.plot(t, T[i,:], label='temperature fluctuation')
plt.xlabel('time (s)')
plt.ylabel('temperature (c)')
plt.legend(loc='lower right')
plt.show()

到目前为止，我已经设法使用热方程的解来使用此代码绘制图形。

import numpy as np
import matplotlib.pyplot as plt
from math import *


t0=28
tm=38
x= 4*10e-2
T=24*3600
w=(2*pi)/T
D=0.75
delta=sqrt(2*D/w)
t=np.arange(0, 232920, 10*60)
temp = t0 + tm*np.exp(-x/delta)*np.cos(w*t-(x/delta))

plt.plot(t, temp, label='temperature fluctuation')
plt.xlabel('time (s)')
plt.ylabel('temperature (c)')
plt.legend(loc='lower right')
plt.show()

我现在正试图让python来解决代码并得到相同的结果