有限元方法求解二维矩形区域椭圆方程

最新推荐文章于 2023-01-12 23:08:49 发布
最新推荐文章于 2023-01-12 23:08:49 发布
阅读量1.7k 收藏 6
点赞数
分类专栏： 传统数值方法 文章标签： 线性代数 算法 矩阵
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/forrestguang/article/details/108428191
版权
传统数值方法 专栏收录该内容
35 篇文章 17 订阅 ¥199.90 ¥99.00
订阅专栏

利用Galerkin方法


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])

class TESET():
    def __init__(self, bounds, nx):
        self.bounds = bounds
        self.nx = nx
        self.hx = [(self.bounds[0,1]-self.bounds[0,0])/self.nx[0],
                   (self.bounds[1,1]-self.bounds[1,0])/self.nx[1]]

        self.size = self.nx[0]*self.nx[1]
        self.X = torch.zeros(self.size,2)
        m = 0
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                self.X[m,0] = self.bounds[0,0] + (i+0.5)*self.hx[0]
                self.X[m,1] = self.bounds[1,0] + (j+0.5)*self.hx[1]
                m = m+1

        self.u_acc = UU(self.X[:,0:2],[0,0]).reshape(self.size,1)
class FENET():
    def __init__(self,bounds,nx):
        self.dim = 2
        self.bounds = bounds
        self.hx = [(bounds[0,1] - bounds[0,0])/nx[0],
                  (bounds[1,1] - bounds[1,0])/nx[1]]
        self.nx = nx
        self.gp_num = 2
        self.gp_pos = [(1 - np.sqrt(3)/3)/2,(1 + np.sqrt(3)/3)/2]
        
        self.Node()
        self.Unit()
        self.matrix()
        
    def Node(self):#生成网格点(M + 1)*(N + 1)
        self.Nodes_size = (self.nx[0] + 1)*(self.nx[1] + 1)
        self.Nodes = torch.zeros(self.Nodes_size,self.dim)
        m = 0
        for i in range(self.nx[0] + 1):
            for j in range(self.nx[1] + 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
    def Unit(self):#生成所有单元，单元数目(M*N)
        self.Units_size = self.nx[0]*self.nx[1]
        self.Units_Nodes = np.zeros([self.Units_size,4],np.int)#每个单元有4个点，记录这4个点的编号
        self.Units_Int_Points = torch.zeros(self.Units_size,#划分成M*N个小区域，每个区域有4个积分点
                                            self.gp_num*self.gp_num,self.dim)
        m = 0
        for i in range(self.nx[0]):
            for j in range(self.nx[1]):
                self.Units_Nodes[m,0] = i*(self.nx[1] + 1) + j
                self.Units_Nodes[m,1] = i*(self.nx[1] + 1) + j + 1
                self.Units_Nodes[m,2] = (i + 1)*(self.nx[1] + 1) + j
                self.Units_Nodes[m,3] = (i + 1)*(self.nx[1] + 1) + j + 1
                n = 0
                for k in range(self.gp_num):
                    for l in range(self.gp_num):
                        self.Units_Int_Points[m,n,0] = \
                            self.bounds[0,0] + (i + self.gp_pos[k])*self.hx[0]
                        self.Units_Int_Points[m,n,1] = \
                            self.bounds[1,0] + (j + self.gp_pos[l])*self.hx[1]
                        n = n + 1
                m = m + 1
    def matrix(self):
        self.A = torch.zeros(self.Nodes_size,self.Nodes_size)#self.Nodes_size = (M+ 1)*(N + 1)
        for m in range(self.Units_size):#self.Units_size = M*N，第m个区域单元的4个积分点
            for k in range(4):
                ind0 = self.Units_Nodes[m,k]#self.Units_Nodes = [M*N,4],第m个区域中第k个网格点，第k个基函数编号
                for l in range(4):
                    ind1 = self.Units_Nodes[m,l]#self.Units_Nodes = [M*N,4],第m个区域中第l网格点，第k个基函数编号
                    #第m个区域上，两个基函数梯度的乘积的积分a(u,v)
                    self.A[ind0,ind1] += self.Int_basic_basic(ind0,ind1,m)
        #边界上的网格点特殊处理
        #左右边界
        for j in range(self.nx[1] + 1):#刚才设置过A的元素，现在需要把对应边界上的元素覆盖掉
            ind = j
            self.A[ind,:] = torch.zeros(1,self.Nodes_size)#self.Nodes_size = (M + 1)*(N + 1)
            self.A[ind,ind] = 1
            ind = self.nx[0]*(self.nx[1] + 1) + j
            self.A[ind,:] = torch.zeros(1,self.Nodes_size)
            self.A[ind,ind] = 1
        #上下边界
        for i in range(1,self.nx[0]):
            ind = i*(self.nx[1] + 1)
            self.A[ind,:] = torch.zeros(1,self.Nodes_size)
            self.A[ind,ind] = 1
            ind = (i + 1)*(self.nx[0] + 1) - 1
            self.A[ind,:] = torch.zeros(1,self.Nodes_size)
            self.A[ind,ind] = 1
    def right(self,UU,FF):
        self.b = torch.zeros(self.Nodes_size,1)#self.Nodes_size = (M + 1)*(N + 1)
        for m in range(self.Units_size):#self.Units_size = M*N
            for k in range(4):
                ind = self.Units_Nodes[m,k]#第m个单元区域内第k个网格点，第k个基函数的编号
                self.b[ind] += self.Int_F_basic(ind,m,FF)
        #边界上的网格点需要特殊处理
        #左右边界
        for j in range(self.nx[1] + 1):# b = [(nx[0] + 1)*(nx[1] + 1),1]
            ind = j#b[i*(nx[1] + 1) + j],b按列存储，这是第一列
            self.b[ind] = UU(self.Nodes[ind:ind + 1,:],[0,0])#self.Nodes存储(M + 1)*(N + 1)个网格点
            ind  = self.nx[0]*(self.nx[1] + 1) + j#这是最后一列元素
            self.b[ind] = UU(self.Nodes[ind:ind + 1,:],[0,0])
        #上下边界
        for i in range(1,self.nx[0]):
            ind = i*(self.nx[1] + 1)#每一列中取第一个元素就是下边界
            self.b[ind] = UU(self.Nodes[ind:ind + 1],[0,0])#self.Nodes存储(M + 1)*(N + 1)个网格点
            ind = (i + 1)*(self.nx[1] + 1) - 1#每一列中取最后一个元素就是上边界
            self.b[ind] = UU(self.Nodes[ind:ind + 1],[0,0])#self.Nodes存储(M + 1)*(N + 1)个网格点
                    
    def phi(self,X,order):#[-1,1]*[-1,1]，在原点取值为1，其他网格点取值为0的基函数
        ind00 = (X[:,0] >= -1);ind01 = (X[:,0] >= 0);ind02 = (X[:,0] >= 1)
        ind10 = (X[:,1] >= -1);ind11 = (X[:,1] >= 0);ind12 = (X[:,1] >= 1)
        if order == [0,0]:
            return (ind00*~ind01*ind10*~ind11).float()*(1 + X[:,0])*(1 + X[:,1]) + \
                    (ind01*~ind02*ind10*~ind11).float()*(1 - X[:,0])*(1 + X[:,1]) + \
                    (ind00*~ind01*ind11*~ind12).float()*(1 + X[:,0])*(1 - X[:,1]) + \
                    (ind01*~ind02*ind11*~ind12).float()*(1 - X[:,0])*(1 - X[:,1])
        if order == [1,0]:
            return (ind00*~ind01*ind10*~ind11).float()*(1 + X[:,1]) + \
                    (ind01*~ind02*ind10*~ind11).float()*(-(1 + X[:,1])) + \
                    (ind00*~ind01*ind11*~ind12).float()*(1 - X[:,1]) + \
                    (ind01*~ind02*ind11*~ind12).float()*(-(1 - X[:,1]))
        if order == [0,1]:
            return (ind00*~ind01*ind10*~ind11).float()*(1 + X[:,0]) + \
                    (ind01*~ind02*ind10*~ind11).float()*(1 - X[:,0]) + \
                    (ind00*~ind01*ind11*~ind12).float()*(-(1 + X[:,0])) + \
                    (ind01*~ind02*ind11*~ind12).float()*(-(1 - X[:,0]))
    def basic(self,X,order,i):#根据网格点的存储顺序，遍历所有网格点，取基函数
        temp = (X - self.Nodes[i,:])/torch.tensor([self.hx[0],self.hx[1]])
        if order == [0,0]:
            return self.phi(temp,order)
        if order == [1,0]:
            return self.phi(temp,order)/self.hx[0]
        if order == [0,1]:
            return self.phi(temp,order)/self.hx[1]
        
    # 计算两个基函数梯度的内积，并积分
    # u_ind：计算积分的单元(同样是索引)
    def Int_basic_basic(self,i,j,u_ind):#表示第i,j个基函数的梯度的乘积，以及在第u_ind个区域的积分
        X0 = self.Units_Int_Points[u_ind,:,:]
        basic0 = torch.zeros(X0.shape)
        basic0[:,0] = self.basic(X0,[1,0],i)
        basic0[:,1] = self.basic(X0,[0,1],i)
        
        X1 = self.Units_Int_Points[u_ind,:,:]
        basic1 = torch.zeros(X1.shape)
        basic1[:,0] = self.basic(X1,[1,0],j)
        basic1[:,1] = self.basic(X1,[0,1],j)
        return ((basic0*basic1).sum(1)).mean()*self.hx[0]*self.hx[1]
    def Int_F_basic(self,i,u_ind,FF):#第i个基函数与右端项乘积，在第u_ind个单元积分
        X = self.Units_Int_Points[u_ind,:,:]
        return (FF(X)*self.basic(X,[0,0],i)).mean()*self.hx[0]*self.hx[1]
    def solve(self,UU,FF):
        self.right(UU,FF)
        self.Nodes_V,lu = torch.solve(self.b,self.A)# 求解方程组得到所有网格点的值
    def Uh(self,X):
        uh = torch.zeros(X.shape[0])
        for i in range(self.Nodes_size):#self.Nodes_size = (M + 1)*(N + 1)
            # 计算数据集 关于 基函数中心 的相对位置
            uh += self.Nodes_V[i]*self.basic(X,[0,0],i)
        return uh.view(-1,1)
    def Ufd(self,X,V):
        uh = torch.zeros(X.shape[0])
        for i in range(V.shape[0]):
            uh += V[i]*self.basic(X,[0,0],i)
        return uh.view(-1,1)
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)
def main():
    st = time.time()
    bounds = torch.tensor([[-1,2],[-1,0]]).float()
    nx_tr = [20,20]    # 网格大小
    nx_te = [100,100]  # 测试集大小
    fe_net = FENET(bounds,nx_tr)
    fe_net.solve(UU,FF)
    
    teset = TESET(bounds,nx_te)
    u_pred = fe_net.Uh(teset.X)
    ERROR = error(u_pred,teset.u_acc).item()
    ela = time.time() - st
    print('the error of FE is %.3e,use time :%.2f'%(ERROR,ela))
if __name__ == '__main__':
    main()
Galerkin码农选手
关注
  • 0
    点赞
  • 踩
  • 6
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 6
    评论
专栏目录
订阅专栏
matlab如何敲代码-EllipticFEM.jl:椭圆，抛物线和双曲型偏微分方程的二维有限元求解器的实现
05-23
matlab如何敲代码EllipticFEM.jl —用Julia编写的椭圆，抛物线和双曲线PDE的FEM解算器 该软件包提供了二维有限元方法的实现，用于求解Julia中的椭圆，抛物线或双曲型偏微分方程。 特别是可以解决 $ (\partial_{t(t)} u) - \nabla (A*\nabla u) + b*\nabla u + c*u = f $ in $\Omega$, $ u = g\_D $ on $\Gamma\_D$, $ (A*\nabla u)*n = g\_N $ on $\Gamma\_N$, $ u = u $ on $\Gamma\_P$, 其中$ \ Omega $是域，而$ \ Gamma_D $，$ \ Gamma_N $和$ \ Gamma_P $分别表示Dirichlet，Neumann和周期边界。 术语$（\ partial_ {t（t）} u）$可以不存在（椭圆方程），也可以是$ \ partial_t u $（抛物线方程）或$ \ partial_ {tt}
用有限元方法求解二维possion方程
12-25
提供了用有限元方法求解三角形定解区域上二维possion方程的matlab程序
二维有限元方程matlab,有限元法求解二维Poisson方程的MATLAB实现
weixin_28850485的博客
03-17 2951
有限元法求解二维 Poisson 方程的 MATLAB 实现 陈 莲a ，郭元辉b ，邹叶童a ( 西华师范大学 a． 数学与信息学院; b． 教育信息技术中心，四川南充 6437009) 摘 要: 文章讨论了圆形区域上的三角形单元剖分、有限元空间，通过变分形式离散得到有限元方程． 用 MATLAB 编程求得数值解，并进行了误差分析． 关键词: Poisson 方程; 有限元方法; MATLAB ...
有限元（FEM) 、有限差分(FDM)和有限体积（FVM) 的优势和劣势
热流体仿真
06-23 3万+
椭圆方程适合用有限元，当然椭圆方程用什么都可以，只不过有限元方法至于椭圆方程就好比青椒对瘦肉，绝配。有限差分百搭，解什么都可以，只不过有限差分大多要求结构网格，一般学术研究中用得多。有限差分还有个不大不小的优势，就是能利用结构网格的拓扑优势轻松搞大模板，构造出高精度格式。有限体积法就比较强大了。除了高精度构造略微麻烦，几乎通吃有限差分所有领域，双曲性、抛物型、椭圆形都可以。当然对于椭圆形方程，不如...
一类二维变系数椭圆方程数值求解
12-30
一类二维变系数椭圆方程数值求解，刘相国，闵涛，利用有限元法求解了二维变系数椭圆方程边值问题，得到了相应的误差分析，并进行了数值模拟。结果表明解此类问题时有限元法具有程
一般情况下的椭圆方程
weixin_46867208的博客
01-12 2833
椭圆的矩阵形式以及一般椭圆方程
poisson（matlab） 一维和二维 有限元程序
11-01
poisson（matlab） 一维和二维 有限元程序，自己初学时写的，相对比较简单，适用于初学者.
有限元矩形网格剖分
04-07
本程序采用MATLAB编写，可以自动进行矩形网格的剖分，以及图形的生成。
有限元方法求解一维扩散方程(FEALPy)
图灵的猫的博客
02-14 3032
有限元方法求解一维扩散方程 之前完成了 FEALPy 有限元求解 Poisson 方程 的数值算例, 通过湘潭大学王唯师兄的协助, 今天基于 FEALPy 运用有限元方法求解一个抛物型方程，为说明使用方法而又简单起见，于是理论部分专注于讨论下面这个一维的抛物型方程
matlab使用有限元方法求解偏微分方程
09-05
大量的有限元法求解偏微分方程的实例程序，程序里有详细的解释语句 大量的有限元法求解偏微分方程的实例程序，程序里有详细的解释语句
四维张量积二次矩形有限元最大模的超逼近 (2011年)
05-20
针对Poisson方程Dirichlet边值问题，首先建立了四维投影型插值算子，并应用它得到了正规剖分下四维张量积二次矩形有限元的弱估计，在此基础上，结合四维离散Green函数的估计，研究四维张量积二次矩形有限元解及梯度...
Matlab求解偏微分方程的代码-SparseGrids:用于通过varionsFEM求解椭圆偏微分方程的MATLAB/Octave程序集
06-01
脚本和函数文件，用于求解单位平方上的线性二维偏微分方程，使用具有双线性元素的标准 Galerkin FEM 和两种稀疏网格方法。 此代码用于在 Russell, S. 和 Madden, N. 中生成结果。稀疏网格有限元方法的分析和实现。 ...
galerkin有限元法matlab实现,有限元法求解二维Poisson方程的MATLAB实现
weixin_35225144的博客
03-18 1034
有限元法求解二维Poisson方程的MATLAB实现陈莲a，郭元辉b，邹叶童a【摘要】文章讨论了圆形区域上的三角形单元剖分、有限元空间，通过变分形式离散得到有限元方程. 用MATLAB编程求得数值解，并进行了误差分析.【期刊名称】洛阳师范学院学报【年(卷),期】2018(037)005【总页数】4【关键词】 Poisson方程；有限元方法； MATLAB编程；三角形单元剖分0 引言热学、流体力学、...
有限元方法入门：有限元方法简单的二维算例（三角形剖分）
热门推荐
LSEC小陆的博客
07-20 4万+
有限元方法简单的二维算例（三角形剖分） 算例描述 我们对下述椭圆边值问题 \label{eq1}{−Δu=fu|∂Ω=0{−Δu=fu|∂Ω=0 \left \{ \begin{aligned} & -\Delta u = f\\ & u|_{\partial \Omega} = 0 & \end{aligned} \right. 其中，Ω=(0,1)2Ω=(0,1)2\Omega = (...
Matlab偏微分方程快速上手：使用pde有限元工具箱求解二维偏微分方程
weixin_47006934的博客
02-01 1万+
Matlab偏微分方程快速上手：使用pdetool工具箱求解二维偏微分方程，适用于数学建模、数学实验，简单的偏微分方程数值计算与工程问题。
二维静电场有限元matlab,实验2 MATLAB有限元计算
weixin_34082139的博客
03-17 1930
MATLAB电磁场有限元计算实验目的：1、了解有限元算法的原理，熟练运用MATLAB环境的PDE工具。 2、熟练运用PDE工具分析简单的电磁场边值问题。实验内容：一、 有限元简介在电磁场的计算中, 仅对那些具有最简单边界条件和场域几何形状规则的问题才有解析解, 多数问题的求解必须用数值计算的方法,其场域分布的数值计算内容是学习难点。本实验将有限元法和Matlab 结合起来对电磁场教学中的电位分布问...
针对neumann边界条件的差分法代码
forrestguang的博客
10-11 4135
需要添加人工的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边界，也就是说矩
五点差分法求解泊松方程-共轭梯度法和G-S迭代法
forrestguang的博客
12-30 2895
差分法求解矩形区域椭圆方程 Δu=f,x∈Ω=[−1,1]∗[−1,1]\Delta u = f, x \in \Omega = [-1,1]*[-1,1]Δu=f,x∈Ω=[−1,1]∗[−1,1] u=g,x∈∂Ωu = g, x \in \partial \Omegau=g,x∈∂Ω 在x轴和y轴上分别取分划长度为dx,dy，令M=int(2/dx)+1,N=int(2/dx)+1M = int(2/dx) + 1,N = int(2/dx) + 1M=int(2/dx)+1,N=int(2/dx)+
matlab编写的求解二维椭圆偏微分方程的一次和二次有限元方法的
最新发布
11-04
二维椭圆偏微分方程是一类常见的数学模型，可以用于描述许多现实世界中的问题，如热传导、流体力学等。而求解这类方程的一次和二次有限元方法是常用的数值求解方法之一。 在Matlab中，可以通过编写相应的程序来实现一次和二次有限元方法的求解过程。首先，需要将二维椭圆偏微分方程离散化为有限元形式，得到一组代数方程。然后，我们可以使用一次或二次有限元方法建立线性或二次插值函数空间，并将离散化后的方程表示为矩阵方程。 对于一次有限元方法，我们使用线性插值函数空间，在Matlab中通常使用`linearshape`函数来定义线性插值。然后，我们可以通过组装刚度矩阵和载荷矩阵得到一个线性方程组，再使用`mldivide`函数求解该方程组。 对于二次有限元方法，我们使用二次插值函数空间，在Matlab中通常使用`quadraticshape`函数来定义二次插值。同样地，我们可以通过组装刚度矩阵和载荷矩阵得到一个二次方程组，再使用`mldivide`函数求解该方程组。 此外，为了更好地进行数值计算，我们还可以使用迭代方法，如共轭梯度法或预处理共轭梯度法来加速求解过程。 总之，通过Matlab编写的一次和二次有限元方法可以较为准确地求解二维椭圆偏微分方程，是一种常用的数值求解方法。在实际应用中，我们还可以通过调节网格密度、选择合适的插值函数和使用更高阶的有限元方法来改进计算结果的精度。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 6
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Galerkin码农选手

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值