为什么解决方案脱离了图形？

Question

晚上好。 我一直在编写一些python代码，使用thomas算法来求解2D-ADI方法的方程，从而求解热方程。 解决该问题的过程如下： 首先，我创建一个获取200x320数据文件的程序（0和1的任意组合都可以，我使用一个提供M形初始轮廓的程序）并为其创建初始条件。 因此，我求解的是200x320点网格（x=200，y=320）

然后，我定义了thomas算法，它解决了所有内部点的方程（它给了我解U，但没有向它添加边界条件，在这种情况下是齐次BCs（=零）。 之后，我对ADI方法本身进行编程，一个从时间级别n到（n+1/2），另一个从（n+1/2）到（n+1），这就是我们想要从n.https://prnt.sc/wjzpp5看到我实现的方程的地方

最后，我想表示9个不同迭代的解决方案，从i=0到1=1500

现在，所有这些都是好的（我认为），因为我得到的解似乎与我对扩散/热方程（缓慢扩散/模糊）的期望一致。我的问题是，当我绘制解矩阵时，解只是…在一些迭代后逃离图片！！！！（https://prnt.sc/wjzsnp，初始轮廓是M）

我的怀疑是，当我存储解决方案时，问题可能在索引中（thomas（a，b，c，d））在每次迭代都将边界条件0添加到矩阵中之后。比如，可能我实际上没有正确地添加边界条件？或者可能我的索引错误，因此，不是每次都将thomas解存储在矩阵的中心，而是将解存储为索引up和索引left-ecery迭代，因此，随着时间的推移，它会“离开”

我正在运行的代码如下，如果有任何建议或任何人发现潜在的问题，我将不胜感激。为我的英语道歉，它不是我的主要语言

import numpy as np
import matplotlib.pyplot as plt
from ThomasF import thomas
from matplotlib import cm

Jx = 200 #Grid w/ 200 points for x
Jy = 320 #Grid w/ 320 points for y
vx = 0.75
vy = 0.75 
#We define the thomas algorithm, which solves for the spatial points except at the boundaries
def thomas(a,b,c,d):
    n = len(b); m = n-1
    C = np.zeros(m); D = np.zeros(n)    

    C[0] = c[0]/b[0]; D[0] = d[0]/b[0]
    for i in range(1, m):
        temp = b[i] - a[i-1]*C[i-1]
        C[i] = c[i]/temp
        D[i] = (d[i] - a[i-1]*D[i-1])/temp
    D[m] = (d[m] - a[m-1]*D[m-1])/(b[m] - a[m-1]*C[m-1])

    x = np.zeros(n); x[m] = D[m]
    for i in range(1, m+1): x[m-i] = D[m-i] - C[m-i]*x[m+1-i]
    return x

#We define the function to read the initial condition file
def llegeixCI():#Reads a 200x320 datapoint file and  creates an initial condition matrix
    file_name="M_200x320.dat"
#Lectura del fitxer línia a línia
    with open(file_name) as f:
        content=f.read().splitlines()
#We define a 200x320 column matrix to store the initial condition
    un=np.zeros((200,320))
#We store the data in the matrix
    for i in range(0,200):
        j=0
    for x in content[i]:
        un[i,j]=x
        j=j+1           
return un
#To start with ( at the beginnint) our profile is the initial condition
U0 = llegeixCI()
#We define the first ADI method, which goes from n to n+1/2, and iterates over x for every y
def ADIx(vx,vy,Jx,Jy,U0):
    a = [*[vx/2] * (Jx-3)]
    b = [*[1+vx] * (Jx-2)]
    c = [*[vx/2] * (Jx-3)]
    d = zeros(Jx-2)
    U = U0
    for s in range(1,Jy-1):
        for r in range(1,Jx-1):
            d[r-1] = 0.5*vy*U[r,s+1]+(1-vy)*U[r,s]+0.5*vy*U[r,s-1]
        U[:,s] = [0,*thomas(a,b,c,d),0] #I store my results in the matrix
    return U
 #We define the second ADI method, which goes from n+1/2 to n+1, and iterates over y for every x
def ADIy(vx,vy,Jx,Jy,U0):
    a = [*[vy/2] * (Jy-3)]
    b = [*[1+vy] * (Jy-2)]
    c = [*[vy/2] * (Jy-3)]
    d = zeros(Jy-2)
    U = ADIx(vx,vy,Jx,Jy,U0)
    for r in range(1,Jx-1):
        for s in range(1,Jy-1):
            d[s-1] = 0.5*vx*U[r+1,s]+(1-vx)*U[r,s]+0.5*vx*U[r-1,s]
        U[r,:] = [0,*thomas(a,b,c,d),0]  #I store my results in the matrix
    return U 

#Graphing
i = 0
I  = iter(range(9))
prints = [0,10,50,100,200,400,700,1000,1500]


fig,ax = plt.subplots(3,3)
ax = ax.flatten()

while i<= 1500:
    if i == 0:
        U0 = llegeixCI()
    else:
         U0 =ADIy(vx,vy,Jx,Jy,U0)
    if i in prints:
        ax[next(I)].imshow(U0,cmap = cm.bone)
        #fig = plt.figure()

    i+=1