针对L型区域的椭圆方程的差分法

#采用五点差分格式，首先给出精确解和右端项的定义。

'''
def UU(X, order):
    temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5
    if order[0]==0 and order[1]==0:
        return torch.log(temp)
    if order[0]==1 and order[1]==0:
        return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
    if order[0]==0 and order[1]==1:
        return temp**(-1) * (20*(X[:,0]+X[:,1]) - 2*(X[:,0]-X[:,1]))
    if order[0]==2 and order[1]==0:
        return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) ** 2 \
               + temp**(-1) * (22)
    if order[0]==0 and order[1]==2:
        return - temp**(-2) * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) ** 2 \
               + temp**(-1) * (22)
'''
def UU(X,order):
    if order == [0,0]:
        return X[:,0]**2 + X[:,1]**2
    if order == [2,0]:
        return 2*torch.ones(X.shape[0])
    if order == [0,2]:
        return 2*torch.ones(X.shape[0])

def FF(X):
    return -UU(X,[2,0]) - UU(X,[0,2])

针对L型区域而言，一共有6个点控制边界，x轴有3个控制点，y轴有3个控制点。给出边界bounds = [[xa,xb,xc],[ya,yb,yc]],例如以下这个L型区域，xa =-1,xb = 0,xc = 1;ya = -1,yb = 0,yc = 1。

##下面重点介绍区域点采样以及矩阵和右端项的构建。

class FENET():
    def __init__(self,bounds,hx):
        self.bounds = bounds
        self.dim = 2
        self.hx = hx
        nx = [[int((bounds[0,1] - bounds[0,0])/hx[0]),
               int((bounds[0,2] - bounds[0,1])/hx[0])],
              [int((bounds[1,2] - bounds[1,0])/hx[1]),
               int((bounds[1,1] - bounds[1,0])/hx[1])]]
        self.nx = torch.tensor(nx)
        
        self.Node()
        self.u_acc = UU(self.Nodes,[0,0]).view(-1,1)
    def Node(self):#生成网格点(M + 1)*(N + 1)
        self.Nodes_size = (self.nx[0,0] + 1)*(self.nx[1,0] + 1) + self.nx[0,1]*(self.nx[1,1] + 1)
        self.Nodes = torch.zeros(self.Nodes_size,self.dim)
        m = 0
        for i in range(self.nx[0,0] + 1):
            for j in range(self.nx[1,0] + 1):
                self.Nodes[m,0] = self.bounds[0,0] + i*self.hx[0]
                self.Nodes[m,1] = self.bounds[1,0] + j*self.hx[1]
                m = m + 1
        for i in range(self.nx[0,1]):
            for j in range(self.nx[1,1] + 1):
                self.Nodes[m,0] = self.bounds[0,1] + (i + 1)*self.hx[0]
                self.Nodes[m,1] = self.bounds[1,0] + j*self.hx[1]
                m = m + 1
    def matrix(self):
        self.A = torch.zeros(self.Nodes_size,self.Nodes_size)
        m = 0
        for i in range(self.nx[0,0] + 1):
            for j in range(self.nx[1,0] + 1):
                dx = self.hx[0];dy = self.hx[1]
                if i == 0 or j == 0 or j == self.nx[1,0]:
                    self.A[m,m] = 1
                    plt.plot(self.Nodes[m,0],self.Nodes[m,1],'r*')
                elif (i == self.nx[0,0] and self.Nodes[m,1] >= self.bounds[1,1]):
                    self.A[m,m] = 1
                    plt.plot(self.Nodes[m,0],self.Nodes[m,1],'r*')
                else:
                    self.A[m,i*(self.nx[1,0] + 1) + j] = 2*(dx/dy + dy/dx)
                    self.A[m,(i - 1)*(self.nx[1,0] + 1) + j] = -dy/dx
                    self.A[m,(i + 1)*(self.nx[1,0] + 1) + j] = -dy/dx
                    self.A[m,i*(self.nx[1,0] + 1) + j - 1] = -dx/dy
                    self.A[m,i*(self.nx[1,0] + 1) + j + 1] = -dx/dy
                m = m + 1
        n = 0
        for i in range(self.nx[0,1]):
            for j in range(self.nx[1,1] + 1):
                dx = self.hx[0];dy = self.hx[1]
                if i == self.nx[0,1] - 1 or j == self.nx[1,1] or j == 0:
                    self.A[m + n,m + n] = 1
                    plt.plot(self.Nodes[m + n,0],self.Nodes[m + n,1],'r*')
                else:
                    self.A[m + n,m + n] = 2*(dx/dy + dy/dx)
                    self.A[m + n,m + (i - 1)*(self.nx[1,1] + 1) + j] = -dy/dx
                    self.A[m + n,m + (i + 1)*(self.nx[1,1] + 1) + j] = -dy/dx
                    self.A[m + n,m + i*(self.nx[1,1] + 1) + j - 1] = -dx/dy
                    self.A[m + n,m + i*(self.nx[1,1] + 1) + j + 1] = -dx/dy
                n = n + 1
    def right(self,UU,FF):
        self.b = torch.zeros(self.Nodes_size,1)
        m = 0
        for i in range(self.nx[0,0] + 1):
            for j in range(self.nx[1,0] + 1):
                dx = self.hx[0];dy = self.hx[1]
                if i == 0 or j == 0 or j == self.nx[1,0]:
                    self.b[m] = UU(self.Nodes[m:m + 1,:],[0,0])
                    plt.plot(self.Nodes[m,0],self.Nodes[m,1],'r*')
                elif (i == self.nx[0,0] and self.Nodes[m,1] >= self.bounds[1,1]):
                    self.b[m] = UU(self.Nodes[m:m + 1,:],[0,0])
                    plt.plot(self.Nodes[m,0],self.Nodes[m,1],'r*')
                else:
                    self.b[m] = dx*dy*FF(self.Nodes[m:m + 1,:])
                m = m + 1
        n = 0
        for i in range(self.nx[0,1]):
            for j in range(self.nx[1,1] + 1):
                dx = self.hx[0];dy = self.hx[1]
                if i == self.nx[0,1] - 1 or j == self.nx[1,1] or j == 0:
                    self.b[m + n] = UU(self.Nodes[m + n:m + n + 1,:],[0,0])
                    plt.plot(self.Nodes[m + n,0],self.Nodes[m + n,1],'r*')
                else:
                    self.b[m + n] = dx*dy*FF(self.Nodes[m + n:m + n + 1,:])
                n = n + 1
    def solve(self,UU,FF):
        self.matrix()
        self.right(UU,FF)
        u_pred,lu = torch.solve(self.b,self.A)
        return u_pred
def error(u_pred,u_acc):
    fenzi = ((u_pred - u_acc)**2).sum()
    fenmu = (u_acc**2 + 1e-8).sum()
    return (fenzi/fenmu)**(0.5)
bounds = torch.tensor([[-2,0,2],[-1,0,1]]).float()
hx = [0.1,0.1]
fenet = FENET(bounds,hx)
u_pred = fenet.solve(UU,FF)
print('the error of FD method on L-square is %.3e'%(error(u_pred,fenet.u_acc)))

这个Node()函数，给出了采样点self.Nodes，以及采样点的数目self.Nodes_size，self.nx[0,0],self.nx[0,1]分别表示x轴上的两个区间[xa,xb],[xb,xc]上的分划数目（注意x轴上的网格点数目为分划数目+1），self.nx[1,0],self.nx[1,1]分别表示x轴上的两个区间[ya,yc],[yc,yb]上的分划数目（注意y轴上的网格点数目为分划数目+2），所以总共有网格点数目(self.nx[0,0] + 1)(self.nx[1,0] + 1) + self.nx[0,1](self.nx[1,1] + 1)。给出采样点的取点顺序以后，后面构建矩阵A的元素就和正常的矩形区域一样，在边界上A[i,i]=1。这里重点就是边界点的索引，为了检查边界点是否取正确，为此我在循环中加入了plt.plot(self.Nodes[m,0],self.Nodes[m,1],‘r*’)，将对应的边界点以图像方式画出来。矩阵A构建成功以后，右端项都是一样的。