针对neumann边界条件的差分法代码

Galerkin码农选手

已于 2022-11-28 19:09:39 修改

阅读量4.6k

点赞数

分类专栏：传统数值方法文章标签： python 深度学习

于 2020-10-11 15:18:34 首次发布

本文链接：https://blog.csdn.net/forrestguang/article/details/109011163

版权

传统数值方法专栏收录该内容

35 篇文章 18 订阅 ¥199.90 ¥299.90

订阅专栏

本文探讨了在处理矩形区域泊松方程时，如何添加人工Dirichlet边界条件以解决纯Neumann边界条件导致的问题。讨论了一阶、二阶差分近似及其局限性，并介绍了构造虚拟点的方法，强调第三种方法效果最佳。同时提到了格心型差分法在处理边界点上的优势。

摘要由CSDN通过智能技术生成

如果不想付费查看博客内容，可以通过本人知乎主页了解添加链接描述

需要添加人工的Dirichlet边界

在这里插入图片描述
考虑的问题仍然是矩形区域的泊松方程。
如果整个区域都采用Neumann边界条件，即 $\partial\Omega_{N}=\partial\Omega$ ，则这个问题没有定解，因为对于某个精确解 $u$ 而言，针对任意的常数C， $u + c$ 都是解
为此，我们常常需要添加一个人工提供的Dirichlet边界条件，比如说在这里，我们选取 $\{-1\}\times[-1,1]$ 为Dirichlet边界，也就是说矩形区域的左边为Dirichlet边界条件。那么问题就修改为：
在这里插入图片描述
假设我们已经对区域进行了对应的网格剖分：
$x_{i}=i*\delta x,i=0,\ldots,M-1$
$y_{j}=j*\delta y,i=0,\ldots,N-1$

下面开始重点阐述在边界上偏导数的处理。

1：一阶差分近似

$\frac{u_{N-1,j} - u_{N-2,j}}{\delta x}=u_{x}(x_{N-1},y_{j})$ ，精度是 $O (h)$ ，其他两条Neumann边界处理一样。

# -*- coding: utf-8 -*-
"""
Created on Sat Oct 10 14:09:35 2020

@author: 2001213226
"""

import torch
import time
import numpy as np
import matplotlib.pyplot as plt
import torch.nn as nn

def UU(X, order,prob):
    if prob==1:
        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]==1 and order[1]==1:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
                   * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
                   + temp**(-1) * (18)
        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)

    if prob==2:
        if order[0]==0 and order[1]==0:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==0:
            return (3*X[:,0]*X[:,0]-1) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==1:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==2 and order[1]==0:
            return (6*X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==1:
            return (3*X[:,0]*X[:,0]-1) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==2:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
    if prob==3:
        temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
        temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
        if order[0]==0 and order[1]==0:
            return temp1 * temp2**(-1)
        if order[0]==1 and order[1]==0:
            return (2*X[:,0]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,0])
        if order[0]==0 and order[1]==1:
            return (-2*X[:,1]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,1])
        if order[0]==2 and order[1]==0:
            return (2) * temp2**(-1) + \
                   2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
        if order[0]==1 and order[1]==1:
            return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
        if order[0]==0 and order[1]==2:
            return (-2) * temp2**(-1) + \
                   2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
def f(prob,X):
    return -UU(X,[2,0],prob)- UU(X,[0,2],prob)
class FD():
    def __init__(self,bound,hx,prob):
        self.prob = prob
        self.dim = 2
        self.hx = hx
        self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,self.dim)
        m = 0
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                self.X[m,0] = bound[0,0] + i*self.hx[0]
                self.X[m,1] = bound[1,0] + j*self.hx[1]
                m = m + 1
        self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
    def matrix(self):
        self.A = torch.zeros(self.nx[0]*self.nx[1],self.nx[0]*self.nx[1])
        dx = self.hx[0];dy = self.hx[1]
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                dx = self.hx[0];dy = self.hx[1]
                if i == 0:
                    
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
                elif i == self.nx[0] - 1:
                    self.A[i*self.nx[1]+j,(i - 1)*self.nx[1]+j] = -1
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
                elif j == 0:
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j + 1] = 1
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = -1
                elif j == self.nx[1] - 1:
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j - 1] = -1
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
                
                elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
                    self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
        return self.A
    def right(self):
        self.b = torch.zeros(self.nx[0]*self.nx[1],1)
        
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                dx = self.hx[0];dy = self.hx[1]
                X = self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:]
                if i == 0:
                    self.b[i*self.nx[1]+j] = UU(X,[0,0],self.prob)
                elif i == self.nx[0] - 1:
                    self.b[i*self.nx[1]+j] = UU(X,[1,0],self.prob)*dx
                elif j == 0:
                    self.b[i*self.nx[1]+j] = UU(X,[0,1],self.prob)*dy
                elif j == self.nx[1] - 1:
                    self.b[i*self.nx[1]+j] = UU(X,[0,1],self.prob)*dy
             
                elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] =  f(self.prob,X)*dx*dy
        return self.b
    def solve(self):
        A = self.matrix()
        b = self.right()
        u,lu = torch.solve(b,A)
        return u
def error(u_pred,u_acc):
    temp = ((u_pred - u_acc)**2).sum()/(u_acc**2).sum()
    return temp**(0.5)
bound = torch.tensor([[0,2],[0,1]]).float()
hx = [0.2,0.1]
prob = 2
fd = FD(bound,hx,prob)
u_pred = fd.solve()
u_acc = fd.u_acc
print(error(u_pred,u_acc))

这个差分法的效果极为不佳。

二阶差分近似

在这里插入图片描述
这里需要注意在 $u_{y}(x,y_{0}),u_{x}(x,y_{N-1})$ 这两个地方做二阶近似的时候，系数相差一个正负号，自己最好推导一遍。

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""

import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn

def UU(X, order,prob):
    if prob==1:
        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]==1 and order[1]==1:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
                   * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
                   + temp**(-1) * (18)
        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)

    if prob==2:
        if order[0]==0 and order[1]==0:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==0:
            return (3*X[:,0]*X[:,0]-1) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==1:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==2 and order[1]==0:
            return (6*X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==1:
            return (3*X[:,0]*X[:,0]-1) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==2:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
    if prob==3:
        temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
        temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
        if order[0]==0 and order[1]==0:
            return temp1 * temp2**(-1)
        if order[0]==1 and order[1]==0:
            return (2*X[:,0]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,0])
        if order[0]==0 and order[1]==1:
            return (-2*X[:,1]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,1])
        if order[0]==2 and order[1]==0:
            return (2) * temp2**(-1) + \
                   2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
        if order[0]==1 and order[1]==1:
            return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
        if order[0]==0 and order[1]==2:
            return (-2) * temp2**(-1) + \
                   2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
def f(prob,X):
    return -UU(X,[2,0],prob)- UU(X,[0,2],prob)
class FD():
    def __init__(self,bound,hx,prob):
        self.prob = prob
        self.dim = 2
        self.hx = hx
        self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,self.dim)
        m = 0
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                self.X[m,0] = bound[0,0] + i*self.hx[0]
                self.X[m,1] = bound[1,0] + j*self.hx[1]
                m = m + 1
        self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
    def matrix(self):
        self.A = torch.zeros(self.size,self.size)
        dx = self.hx[0];dy = self.hx[1]
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                if i == 0:
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
                    
                elif i == self.nx[0] - 1:
                    self.A[i*self.nx[1]+j,(i - 2)*self.nx[1]+j] = 1
                    self.A[i*self.nx[1]+j,(i - 1)*self.nx[1]+j] = -4
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 3
                elif j == 0:
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j + 2] = -1
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j + 1] = 4
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = -3
                elif j == self.nx[1] - 1:
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j - 2] = 1
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j - 1] = -4
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 3
                elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
                    self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
        return self.A
    def right(self):
        self.b = torch.zeros(self.size,1)
        dx = self.hx[0];dy = self.hx[1]
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                X = self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:]
                if i == 0:
                    self.b[i*self.nx[1]+j] = UU(X,[0,0],self.prob)
                    #self.b[i*self.nx[1]+j] = UU(X,[1,0],self.prob)*2*dx 
                elif i == self.nx[0] - 1:
                    self.b[i*self.nx[1]+j] = UU(X,[1,0],self.prob)*2*dx
                elif j == 0:
                    self.b[i*self.nx[1]+j] = UU(X,[0,1],self.prob)*2*dy
                elif j == self.nx[1] - 1:
                    self.b[i*self.nx[1]+j] = UU(X,[0,1],self.prob)*2*dy
                elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] =  f(self.prob,X)*dx*dy
        return self.b
    def solve(self):
        A = self.matrix()
        b = self.right()
        u,lu = torch.solve(b,A)
        return u
def error(u_pred,u_acc):
    temp = ((u_pred - u_acc)**2).sum()/(u_acc**2).sum()
    return temp**(0.5)
bound = torch.tensor([[0,2],[0,1]]).float()
hx = [0.2,0.1]
prob = 3
fd = FD(bound,hx,prob)
u_pred = fd.solve()
u_acc = fd.u_acc
print(error(u_pred,u_acc))

3：构造虚拟点

在这里插入图片描述
这种近似方法的误差也是 $O(h^{2})$ ，但是构造特别麻烦，在矩形区域的4个顶点上需要和其他Neumann边界进行重新划分



import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn

def UU(X, order,prob):
    if prob==1:
        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]==1 and order[1]==1:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
                   * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
                   + temp**(-1) * (18)
        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)

    if prob==2:
        if order[0]==0 and order[1]==0:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==0:
            return (3*X[:,0]*X[:,0]-1) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==1:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==2 and order[1]==0:
            return (6*X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==1:
            return (3*X[:,0]*X[:,0]-1) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==2:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
    if prob==3:
        temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
        temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
        if order[0]==0 and order[1]==0:
            return temp1 * temp2**(-1)
        if order[0]==1 and order[1]==0:
            return (2*X[:,0]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,0])
        if order[0]==0 and order[1]==1:
            return (-2*X[:,1]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,1])
        if order[0]==2 and order[1]==0:
            return (2) * temp2**(-1) + \
                   2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
        if order[0]==1 and order[1]==1:
            return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
        if order[0]==0 and order[1]==2:
            return (-2) * temp2**(-1) + \
                   2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
def f(prob,X):
    return -UU(X,[2,0],prob)- UU(X,[0,2],prob)
class FD():
    def __init__(self,bound,hx,prob):
        self.prob = prob
        self.dim = 2
        self.hx = hx
        self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,self.dim)
        m = 0
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                self.X[m,0] = bound[0,0] + i*self.hx[0]
                self.X[m,1] = bound[1,0] + j*self.hx[1]
                m = m + 1
        self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
    def matrix(self):
        self.A = torch.zeros(self.size,self.size)
        dx = self.hx[0];dy = self.hx[1]
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                if i == 0:
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
                
                elif (i == self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
                    
                elif (j == 0 and i > 0 and i < self.nx[0] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -2*dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
                    self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
                    
                elif (j == self.nx[1] - 1 and i > 0 and i < self.nx[0] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -2*dx/dy
                    
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
                    self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
                
                elif (i == self.nx[0] - 1 and j == 0):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -2*dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
                elif (i == self.nx[0] - 1 and j == self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -2*dx/dy
                    
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
                elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
                    self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
                    
                
        return self.A
    def right(self):
        self.b = torch.zeros(self.size,1)
        dx = self.hx[0];dy = self.hx[1]
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                X = self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:]
                if i == 0:
                    self.b[i*self.nx[1]+j] = UU(X,[0,0],self.prob)
                elif (i == self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy
                elif (j == 0 and i > 0 and i < self.nx[0] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy - UU(X,[0,1],self.prob)*2*dx
                elif (j == self.nx[1] - 1 and i > 0 and i < self.nx[0] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[0,1],self.prob)*2*dx
                
                
                elif (i == self.nx[0] - 1 and j == 0):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy - UU(X,[0,1],self.prob)*2*dx
                elif (i == self.nx[0] - 1 and j == self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy + UU(X,[0,1],self.prob)*2*dx
                elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] =  f(self.prob,self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:])*dx*dy
        return self.b
    def solve(self):
        A = self.matrix()
        b = self.right()
        u,lu = torch.solve(b,A)
        return u
def error(u_pred,u_acc):
    temp = ((u_pred - u_acc)**2).sum()/(u_acc**2).sum()
    return temp**(0.5)
bound = torch.tensor([[0,2],[0,1]]).float()
hx = [0.2,0.1]
prob = 2
fd = FD(bound,hx,prob)
u_pred = fd.solve()
u_acc = fd.u_acc
print(error(u_pred,u_acc))

综合来看，第三种方法的效果最好。

格心型差分法

刚才提到的3种都是格点型差分法代码，格点型指的是取网格单元的顶点来计算，格心型差分法代码指的是取网格中心点来计算，比如说参考下面这个图
在这里插入图片描述
红色的点就是我们采用格点型差分法需要采集的样本点，绿色的点就是格心型我们需要采集的样本点。对于格心型差分法，在内点的处理与格点型没有区别（这里格心型的内点指的是这些网格点的内部，比如说图中的内点只有两个，其余7个都算外点），对于外点也就是边界点，格心型的处理方法就是上述第3种处理，这样最方便，事实上，好像也没得选。下面看代码

import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn
def UU(X, order,prob):
    if prob==1:
        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]==1 and order[1]==1:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
                   * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
                   + temp**(-1) * (18)
        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)

    if prob==2:
        if order[0]==0 and order[1]==0:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==0:
            return (3*X[:,0]*X[:,0]-1) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==1:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==2 and order[1]==0:
            return (6*X[:,0]) * \
                   0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==1:
            return (3*X[:,0]*X[:,0]-1) * \
                   (torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==2:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
    if prob==3:
        temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
        temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
        if order[0]==0 and order[1]==0:
            return temp1 * temp2**(-1)
        if order[0]==1 and order[1]==0:
            return (2*X[:,0]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,0])
        if order[0]==0 and order[1]==1:
            return (-2*X[:,1]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,1])
        if order[0]==2 and order[1]==0:
            return (2) * temp2**(-1) + \
                   2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
        if order[0]==1 and order[1]==1:
            return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
        if order[0]==0 and order[1]==2:
            return (-2) * temp2**(-1) + \
                   2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
def f(prob,X):
    return -UU(X,[2,0],prob)- UU(X,[0,2],prob)
class FD():
    def __init__(self,bound,hx,prob):
        self.prob = prob
        self.dim = 2
        self.hx = hx
        self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]),int((bound[1,1] - bound[1,0])/self.hx[1])]
        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,self.dim)
        m = 0
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                self.X[m,0] = bound[0,0] + (i + 0.5)*self.hx[0]
                self.X[m,1] = bound[1,0] + (j + 0.5)*self.hx[1]
                m = m + 1
        self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
    def matrix(self):
        self.A = torch.zeros(self.size,self.size)
        dx = self.hx[0];dy = self.hx[1]
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                if i == 0 and j == 0:
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
                elif (i == 0 and j > 0 and j < self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
                    self.A[i*self.nx[1]+j,(i + 1)*self.nx[1]+j] = -2*dy/dx
                elif (i == 0 and j == self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -2*dx/dy
                    
                    self.A[i*self.nx[1]+j,(i + 1)*self.nx[1]+j] = -2*dy/dx
                elif (i == self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
                    
                elif (j == 0 and i > 0 and i < self.nx[0] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -2*dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
                    self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
                    
                elif (j == self.nx[1] - 1 and i > 0 and i < self.nx[0] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -2*dx/dy
                    
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
                    self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
                
                elif (i == self.nx[0] - 1 and j == 0):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -2*dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
                elif (i == self.nx[0] - 1 and j == self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -2*dx/dy
                    
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
                elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
                    self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
                    self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
                    self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
                    
                
        return self.A
    def right(self):
        self.b = torch.zeros(self.size,1)
        dx = self.hx[0];dy = self.hx[1]
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                X = self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:]
                if i == 0 and j == 0:
                    self.b[i*self.nx[1]+j] = UU(X,[0,0],self.prob)
                elif (i == 0 and j > 0 and j < self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy - UU(X,[1,0],self.prob)*2*dy
                elif (i == 0 and j == self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy - UU(X,[1,0],self.prob)*2*dy + UU(X,[0,1],self.prob)*2*dx
                elif (i == self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy
                elif (j == 0 and i > 0 and i < self.nx[0] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy - UU(X,[0,1],self.prob)*2*dx
                elif (j == self.nx[1] - 1 and i > 0 and i < self.nx[0] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[0,1],self.prob)*2*dx
                
                
                elif (i == self.nx[0] - 1 and j == 0):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy - UU(X,[0,1],self.prob)*2*dx
                elif (i == self.nx[0] - 1 and j == self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy + UU(X,[0,1],self.prob)*2*dx
                elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
                    self.b[i*self.nx[1]+j] =  f(self.prob,self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:])*dx*dy
        return self.b
    def solve(self):
        A = self.matrix()
        b = self.right()
        u,lu = torch.solve(b,A)
        return u
def error(u_pred,u_acc):
    temp = ((u_pred - u_acc)**2).sum()/(u_acc**2).sum()
    return temp**(0.5)
def LE(u_pred,u_acc):
    temp = abs(u_pred - u_acc)
    return max(temp)
bound = torch.tensor([[0,1],[0,1]]).float()
hx = [0.1,0.1]
prob = 3
fd = FD(bound,hx,prob)
u_pred = fd.solve()
u_acc = fd.u_acc
print('prob = %d,error = %.4f'%(prob,error(u_pred,u_acc)))
print(LE(u_pred,u_acc))

Galerkin码农选手

关注

0
点赞
踩
12

收藏

觉得还不错? 一键收藏
打赏
3
评论
针对neumann边界条件的差分法代码

需要添加人工的Dirichlet边界考虑的问题仍然是矩形区域的泊松方程。如果整个区域都采用Neumann边界条件，即∂ΩN=∂Ω\partial\Omega_{N}=\partial\Omega∂ΩN=∂Ω，则这个问题没有定解，因为对于某个精确解uuu而言，针对任意的常数C，u+cu+cu+c都是解为此，我们常常需要添加一个人工提供的Dirichlet边界条件，比如说在这里，我们选取{−1}×[−1,1]\{-1\}\times[-1,1]{−1}×[−1,1]为Dirichlet边界，也就是说矩
复制链接

扫一扫