为什么解决方案会逃离图表？

编程问答 2022-05-31

问题描述

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

然后，我定义了 thomas 算法，该算法解决了所有内部点的方程（它给了我解 U 但不向它添加边界条件，在这种情况下是齐次 BC（= 零）。 之后，我自己编写 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）

我在这里怀疑问题可能出在索引中，当我将边界条件 0 添加到矩阵后，为它经历的每次迭代存储解决方案 (thomas(a,b,c,d))。就像，也许我实际上没有正确添加边界条件？或者也许我的索引错误，所以不是每次都将 thomas 解决方案存储在矩阵的中心，而是将解决方案存储在一个索引向上和左侧 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,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,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,U0):
    a = [*[vy/2] * (Jy-3)]
    b = [*[1+vy] * (Jy-2)]
    c = [*[vy/2] * (Jy-3)]
    d = zeros(Jy-2)
    U = ADIx(vx,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,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,U0)
    if i in prints:
        ax[next(I)].imshow(U0,cmap = cm.bone)
        #fig = plt.figure()

    i+=1

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

indexingloopsnumerical-methodsphysicspython