使用隐式方法Python 1d PDE

Question

我如何编码以下问题的方程以及初始和边界条件？我手工做了它们，但我不确定如何编码它们。我还附加了到目前为止的代码。我知道我在下面编码的内容是不正确的，但是有些相似。

enter Image Description

import numpy as np
import matplotlib.pyplot as plt

# Constants
r1 = 0.06 # Inner radius (m)
r2 = 0.15 # Outer radius (m)
r = r2-r1
a1 = 45.4*10**-6 # Inner thermal diffusivity (m^2/s)
a2 = 101.2*10**-6 # Outer thermal diffusivity (m^2/s)
a = (a1+a2)/2
T1 = 303 # Inner uniform temperature (K), where t = 0 s
T2 = 393 # Outer uniform temperature (K), where t = 0 s
tt = 200 # Total time (s)
dr = 0.001 # (m)
dt = 0.05 # (s)
tol = 0.1 # Absolute tolerance (K)
lam = a*dt/dr/dr # Parameter for convergence rate
print('lambda = % f ' %(lam));
nt = int(tt/dt)+1 # Number of time segments
nr = int(r/dr)+1
rr = np.linspace(0,r,num=nr,endpoint=True)
nn = nr-1 # One fixed boundary
# Create tridiagonal matrix
Ea = np.ones(nn)
Eb = np.ones(nn-1)
# Implicit euler
IE = np.diagflat(-lam*Eb,k=-1) \
+ np.diagflat((1+2*lam)*Ea,k=0) \
+ np.diagflat(-lam*Eb,k=1);
IE[-1,-2] = -2*lam;
# Vector of constants via BCs
ie = np.zeros(nn);
ie[0] = -lam*T1;
# Space and time matrix for T
M = T1*np.ones(nn); # Initial condition
# Implicit
MM_ie = np.zeros((nn,nt))
MM_ie[:,0] = M
# Solve inverse mats here
IE_inv = np.linalg.inv(IE)
for i in range(0,nt-1):
    MM_ie[:,i+1] = np.matmul(IE_inv,MM_ie[:,i]-ie);
wa = T1*np.ones((1,nt));
MM_ie = np.concatenate((wa,MM_ie))

def show_plot(time) :
    tind = int(time/dt)
    print("time index: %d" % tind)
    print(MM_ie.shape)
    fig = plt.figure()
    ax = plt.axes(xlim=(0, r), ylim=(0, 2*T1))
    ax.plot(rr, MM_ie[:,tind], lw=2)
    time_template = 'time = %4.2f s'
    time_text = ax.text(0.75, 0.90,time_template % time,transform=ax.transAxes)
    plt.xlabel('Radius (m)')
    plt.ylabel('Temperature, K')
    plt.title('Temperature Profile')
    plt.show()
show_plot(1)
show_plot(10)
show_plot(30)
show_plot(70)
show_plot(200)