精华内容
下载资源
问答
  • 求解抛物型方程的紧差分方程的加速并行迭代法,郭瑜超,,本文主要研究了抛物方程的紧式差分方程的加速分段并行迭代法。证明了迭代的收敛性,推出了网格加密时迭代的收敛性质。理论分析表
  • 求解双曲方程的隐式差分方程的加速并行迭代法,郭瑜超,,本文主要研究双曲方程隐式差分格式并行迭代的加速方法,算法改进了分裂矩阵的方法,提高了迭代的收敛速度。理论分析证明了它的收
  • 有限差分法求解偏微分方程

    万次阅读 多人点赞 2016-11-06 14:23:03
    只含有未知多元函数及其偏导数的方程,称之为偏微分方程。偏微分方程在工程技术(甚至图像处理)中有...本文主要介绍利用有限差分(涉及Jacobi迭代和Gauss-Siedel迭代求解偏微分方程(以Laplace方程为例)的方法,

    自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元函数及其偏导数的方程,称之为偏微分方程。初始条件和边界条件称为定解条件,未附加定解条件的偏微分方程称为泛定方程。对于一个具体的问题,定解条件与泛定方程总是同时提出。定解条件与泛定方程作为一个整体,称为定解问题。

    在介绍有限差分法求解偏微分方程的过程中,我们会用到Jacobi迭代法与Gauss-Seidel迭代法的相关内容,如果读者对此不甚了解,可以参阅下文:

    欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji,为保证公式、图表得以正确显示,强烈建议你从该地址上查看原版博文。本博客主要关注方向包括:数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。


    椭圆方程

    由于大多数工程问题都是二维问题,所以得到的微分方程一般都是偏微分方程,对于一维问题得到的是常微分方程,解法与偏微分方程类似,故为了不是一般性,这里只讨论偏微分方程。由于工程中高阶偏微分较少出现,所以本文仅仅给出二阶偏微分方程的一般形式,对于高阶的偏微分,可进行类似地推广。二阶偏微分方程的一般形式如下:

    AΦxx+BΦxy+CΦyy=f(x,y,Φ,Φx,Φy)

    其中,Φ 表示一个连续函数。当A,B,C都是常数时,上式成为准线性,有三种准线性方程形式:

    • 如果Δ=B24AC<0,则称为椭圆型方程;
    • 如果Δ=B24AC=0,则称为抛物型方程;
    • 如果Δ=B24AC>0,则称为双曲型方程;

    椭圆方程是工程技术应用中所涉及的偏微分方程里最为普遍的一种形式。根据椭圆方程的具体形式又可以将其分为以下三种形式:

    • 拉普拉斯(Laplace)方程:2u=0
    • 泊松(Poisson)方程:2u=g(x,y)
    • 亥姆霍兹(Helmholtz)方程:2u+λu=0

    其中u是关于xy的二元函数。


    有限差分法

    差分方法又称为有限差分方法或网格法,是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是:先对求解区域作网格剖分,将自变量的连续变化区域用有限离散点(网格点)集代替;将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替;通过用网格点上函数的差商代替导数,将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组(称为差分格式)。如果差分格式有解,且当网格无限变小时其解收敛于原微分方程定解问题的解,则差分格式的解就作为原问题的近似解(数值解)。因此,用差分方法求偏微分方程定解问题一般需要解决以下问题:

    1. 选取网格;
    2. 对微分方程及定解条件选择差分近似,列出差分格式;
    3. 求解差分格式;
    4. 讨论差分格式解对于微分方程解的收敛性及误差估计。

    下面我们就以拉普拉斯方程的数值解法为例来演示一下有限差分法的基本思路。我们首先写出完整的拉普拉斯方程如下:

    2fx2+2fy2=0

    现在的问题其实是要求我们在一个给定的二维区域中求解满足方程的每一点(x,y)。一些区域中的点将被用来给出边界条件(hold boundary conditions)。

    于是我们将整个二维区域离散化成若干个点,如下图所示为其中的五个相邻点:



    根据偏导数的定义则有:
    fx1Δ[f(x+Δ2,y)f(xΔ2,y)]

    2fx21Δ{1Δ[f(x+Δ,y)f(x,y)]1Δ[f(x,y)f(xΔ,y)]}=1Δ2[f(x+Δ,y)2f(x,y)+f(xΔ,y)]

    同理可得:
    2fy21Δ2[f(x,y+Δ)2f(x,y)+f(x,yΔ)]

    将上述两个结果带入拉普拉斯方程可得:
    1Δ2[f(x+Δ,y)+f(xΔ,y)+f(x,y+Δ)+f(x,yΔ)4f(x,y)]=0


    方程组求解

    回想雅各比迭代法(可以参考本文最开始给出的文章链接),假设我们有一个由n个线性方程组成的系统(也就是线性方程组):

    Ax=b

    那么Jacobi迭代可以描述为:
    xki=1ai,i[bijiai,jxk1j]

    其中上标k表示第k轮迭代。

    注意我们在上一小节最后得出的拉普拉斯方程离散化形式给出了(离散化后)区域上众多点中的一个点的求解方程,所有点的求解方程合在一起就构成了一个大的方程组。我们把求解某点(x,y)的方程重写成Jacobi迭代的形式,则有:

    fk(x,y)=14[fk1(xΔ,y)+fk1(x,yΔ)+fk1(x+Δ,y)+fk1(x,y+Δ)]

    重复应用上述迭代式,最后方程就会收敛到解的附近。
    本来连续的一个区域经过离散化处理之后就变成了一个网格结构,假设网格的大小是n2,它们的标签为x1,x2,,xn2,如下图所示



    上面这种自然排列的点序可以得出不超过n2个五元线性方程:
    xin+xi14xi+xi+1+xi+n=0

    注意我们不对处于边界上的点(例如:x1,x2,,x4,x5,x8,x9,)应用上述方程。最后我们将得到如下所示的一个大型的稀疏线性方程组。


    除了Jacobi迭代之外,(如果你看了本文最开始推荐的文章也应该知道)我们还可以采用Guass-Siedel迭代来加速方程组解的收敛速度。Guass-Siedel迭代的一般形式为:
    xki=1ai,i[bij=1i1ai,jxkjj=i+1nai,jxk1j]

    此时拉普拉斯方程需用下式进行求解(我们不再做详细讨论):
    fk(x,y)=14[fk(xΔ,y)+fk(x,yΔ)+fk1(x+Δ,y)+fk1(x,y+Δ)]


    一个讨论:这些内容有什么用?

    我们前面提到偏微分方程在工程中有重要应用。但是在信息技术中有没有离我们比较近的应用实例呢?事实上,偏微分方程在图像处理中就有重要应用!本文主要是以拉普拉斯方程的数值解为例来讨论的,而本文前面我们也提到过椭圆方程中除了拉普拉斯方程之外,还有一类叫做泊松方程。图像处理中基于泊松方程的算法构成了一大类的具有广泛应用的算法,可以用于图像融合、图像去雾、图像拼接等等。例如下图就是基于解泊松方程的方法实现的图像泊松编辑的效果图:


    更多关于泊松融合、泊松编辑方面的内容还可以参考如下链接:
    http://blog.csdn.net/baimafujinji/article/details/50583497

    更多关于数学在图像处理方面的应用,或者图像处理中所需的数学知识和数学原理,你还可以参考我的新书《图像处理中的数学修炼》



    参考文献与推荐阅读材料

    【1】张文生. 科学计算中的偏微分方程有限差分法,高等教育出版社,2006.
    【2】左飞. 图像处理中的数学修炼,清华大学出版社,2016.
    【3】百度文库:有限差分法求解偏微分方程,最后浏览时间为2016年11月.

    展开全文
  • 差分法求解矩形区域椭圆方程 Δ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轴上分别取分划...

    差分法求解矩形区域椭圆方程

    Δu=f,xΩ=[1,1][1,1]\Delta u = f, x \in \Omega = [-1,1]*[-1,1]
    u=g,xΩu = g, x \in \partial \Omega

    在x轴和y轴上分别取分划长度为dx,dy,令M=int(2/dx)+1,N=int(2/dx)+1M = int(2/dx) + 1,N = int(2/dx) + 1,则得到(M + 1)(N + 1)个网格点。假设u在网格点(x_i,y_j)的近似解为u_{i,j}。则对于内部网格点,有:
    (ui+1,j2ui,j+ui1,j)/(dx2)(ui,j+12ui,j+ui,j1)/(dy2)=fi,j-(u_{i+1,j} - 2u_{i,j} + u_{i -1,j})/(dx**2)-(u_{i,j+1} - 2u_{i,j} + u_{i ,j-1})/(dy**2) = f_{i,j}
    两边同时乘dx
    dy,则有:
    (dy/dx)(ui+1,j+ui1,j)2(dx/dy+dy/dx)ui,j+(dx/dy)(ui,j+1+ui,j1)=fi,jdxdy(dy/dx)*(u_{i+1,j}+u_{i-1,j}) -2(dx/dy+dy/dx)u_{i,j}+(dx/dy)*(u_{i,j+1}+u_{i,j-1}) = f_{i,j}*dx*dy

    第一种:生成矩阵求解

    def UU(X, order,prob):#X表示(x,t)
        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:#对x求偏导
                return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
            if order[0]==0 and order[1]==1:#对t求偏导
                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 FF(prob,X):
        return -UU(X,[0,2],prob) - UU(X,[2,0],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],self.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 or i == self.nx[0] - 1 or j == 0 or j == self.nx[1] - 1:
                        self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
                    
                    else:
                        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 = i*dx;y = j*dy
                    X = torch.tensor([[x,y]]).float()
                    if i== 0 or i == self.nx[0] - 1 or j == 0 or j == self.nx[1] - 1:
                        self.b[i*self.nx[1]+j] = UU(self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:],[0,0],self.prob)
                 
                    else:
                        self.b[i*self.nx[1]+j] =  FF(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().numpy()
            b = self.right().numpy()
            x = np.zeros([b.shape[0],1])
            u = GS(A,b,x,200)
            return u
    def error(u_pred,u_acc):
        temp = ((u_pred - u_acc.numpy())**2).sum()/(u_acc.numpy()**2).sum()
        return temp**(0.5)
    bound = torch.tensor([[0,2],[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(error(u_pred,u_acc))
    
    
    

    上面这个方法最原始,同时最耗时,尤其是求解大规模矩阵的时候,矩阵的存储求解特别困难,上面是运用高斯迭代求解。基于这个,我们开始引入线高斯迭代和共轭梯度法,其中最重要的是共轭梯度法求解,这个速度非常惊人。

    线高斯

    import numpy as np
    import time
    def TH(d,l,u,b):
        n = b.shape[0]
        y = np.zeros([n,1]);x = np.zeros([n,1])
        alpha = np.zeros_like(d)
        beta = np.zeros_like(u)
        gama = l.copy()
        alpha[0,0] = d[0,0];beta[0,0] = u[0,0]/d[0,0]
        for i in range(1,n - 1):
            alpha[i,0] = d[i,0] - l[i - 1,0]*beta[i - 1,0]
            beta[i,0] = u[i,0]/alpha[i,0]
        alpha[n - 1,0] = d[n - 1,0] - l[n - 2,0]*beta[n - 2,0]
        y[0,0] = b[0,0]/alpha[0,0]
        for i in range(1,n):
            y[i,0] = (b[i,0] - gama[i - 1,0]*y[i - 1,0])/alpha[i,0]
        x[n - 1,0] = y[n - 1,0]
        for j in range(n - 2,-1,-1):
            x[j,0] = y[j,0] - beta[j,0]*x[j + 1,0]
        return x
    
    def UU(X, order,prob):#X表示(x,t)
        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 np.log(temp)
            if order[0]==1 and order[1]==0:#对x求偏导
                return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
            if order[0]==0 and order[1]==1:#对t求偏导
                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*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
            if order[0]==1 and order[1]==0:
                return (3*X[:,0]*X[:,0]-1) * \
                       0.5*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
            if order[0]==0 and order[1]==1:
                return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                       (np.exp(2*X[:,1])-np.exp(-2*X[:,1]))
            if order[0]==2 and order[1]==0:
                return (6*X[:,0]) * \
                       0.5*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
            if order[0]==1 and order[1]==1:
                return (3*X[:,0]*X[:,0]-1) * \
                       (np.exp(2*X[:,1])-np.exp(-2*X[:,1]))
            if order[0]==0 and order[1]==2:
                return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                       2*(np.exp(2*X[:,1])+np.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 FF(prob,X):
        return -UU(X,[0,2],prob) - UU(X,[2,0],prob)
    np.random.seed(1234)
    
    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 = np.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
        def u_init(self):
            u = np.zeros([self.nx[0],self.nx[1]])
            x = self.X.reshape([self.nx[0],self.nx[1],self.dim])
            u[0,:] = UU(x[0,:,:],[0,0],self.prob)
            u[-1,:] = UU(x[-1,:,:],[0,0],self.prob)
            u[1:self.nx[0] - 1,0] = UU(x[1:self.nx[0] - 1,0,:],[0,0],self.prob)
            u[1:self.nx[0] - 1,-1] = UU(x[1:self.nx[0] - 1,-1,:],[0,0],self.prob)
            return u
        def solve(self,rig):
            tic = time.time()
            u_old = self.u_init()
            u_new = u_old.copy()
            M = self.nx[0];N = self.nx[1]
            dx = self.hx[0];dy = self.hx[1]
            right = (dx*dy*rig).reshape([M - 2,N - 2])
            #print(right[0,:].shape,u_new[0,1:N - 2].shape)
            right[0,:] += u_new[0,1:N - 1]*dy/dx
            
            right[-1,:] += u_new[- 1,1:N - 1]*dy/dx
            
            r1 = dy/dx;r2 = dx/dy
            l = - np.ones([M - 3,1])*r1
            u = - np.ones([M - 3,1])*r1
            d = np.ones([M - 2,1])*2*(r1 + r2)
            #print(d.shape)
            for k in range(1000):
                for j in range(1,N - 1):
                    #print(u_old[1:M - 1,j - 1].shape,right[:,j].shape)
                    b = r2*(u_new[1:M - 1,j - 1] + u_old[1:M - 1,j + 1]) + right[:,j - 1]
                    
                    u_new[1:M - 1,j:j + 1] = TH(d,l,u,b.reshape(-1,1))
                #print(np.linalg.norm(u_new - u_old))
                if (np.linalg.norm(u_new - u_old) < 1e-7):
                    break
                else:
                    u_old = u_new.copy()
                if k%100 == 0:
                    print('the iteration = %d'%(k + 1))
            ela = time.time() - tic
            print('the end iteration is %d,the time:%.2f'%(k + 1,ela))
            return u_new.reshape(-1,1)
    bound = np.array([[0,2.0],[0,1.0]])
    hx = [0.05,0.05]
    prob = 1
    fd = FD(bound,hx,prob)
    M = fd.nx[0];N = fd.nx[1]
    X = fd.X
    u_acc = UU(X,[0,0],prob).reshape(-1,1)
    rig_in = (FF(prob,X).reshape(M,N))[1:M - 1,1:N - 1]
    u_pred = fd.solve(rig_in)
    
    print(max(abs(u_acc - u_pred)))
    
    

    上面这个线高斯迭代避免了矩阵的生成,可以有效加快速度,但是还是太慢了,不能求解大规模问题。下面重点看共轭梯度法,这个的求解速度惊人,可以有效求解大规模问题,同时避免系数矩阵的存储和参与运算。

    共轭梯度法

    import torch
    import numpy as np
    import time
    def UU(X, order,prob):#X表示(x,t)
        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 np.log(temp)
            if order[0]==1 and order[1]==0:#对x求偏导
                return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
            if order[0]==0 and order[1]==1:#对t求偏导
                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)
    def FF(prob,X):
        return -UU(X,[0,2],prob) - UU(X,[2,0],prob)
    
    
    class MG():
        def __init__(self,bound,prob):
            self.prob = prob
            self.dim = 2
            self.bound = bound
        def u_init(self,hx):
            nx = [int((self.bound[0,1] - self.bound[0,0])/hx[0]) + 1,int((self.bound[1,1] - self.bound[1,0])/hx[1]) + 1]
            size = nx[0]*nx[1]
            X = np.zeros([size,self.dim])
            m = 0
            for i in range(nx[0]):
                for j in range(nx[1]):
                    X[m,0] = self.bound[0,0] + i*hx[0]
                    X[m,1] = self.bound[1,0] + j*hx[1]
                    m = m + 1
            x = X.reshape([nx[0],nx[1],self.dim])
            u = np.zeros([nx[0],nx[1]])
            u[0,:] = UU(x[0,:,:],[0,0],self.prob)
            u[-1,:] = UU(x[-1,:,:],[0,0],self.prob)
            u[1:nx[0] - 1,0] = UU(x[1:nx[0] - 1,0,:],[0,0],self.prob)
            u[1:nx[0] - 1,-1] = UU(x[1:nx[0] - 1,-1,:],[0,0],self.prob)
            return u
        def residual(self,u,rig,epoch):#rig也是一个矩阵,跟hx有关,hx控制rig的大小,扁平化处理以后就是方程的右端项
            #由于u已经初始化,把边界条件吸收进去,所以rig的边界为0
            tic = time.time()
            nx = [rig.shape[0],rig.shape[1]]
            hx = [(self.bound[0,1] - self.bound[0,0])/(nx[0] - 1),(self.bound[1,1] - self.bound[1,0])/(nx[1] - 1)]
            
            M = nx[0];N = nx[1]
            r1 = hx[1]/hx[0];r2 = hx[0]/hx[1]
            rig = rig*hx[0]*hx[1]
    
            #共轭梯度法求解线性方程组
            res = np.zeros_like(rig)
            for j in range(1,N - 1):
                res[1:M - 1,j] = rig[1:M - 1,j] - 2*(r1 + r2)*u[1:M - 1,j] + \
                r2*(u[1:M - 1,j - 1] + u[1:M - 1,j + 1]) + \
                r1*(u[0:M - 2,j] + u[2:M,j])
            rho = (res*res).sum()
            
            k = 0
            p = np.zeros_like(res)
            while np.sqrt(rho) > 1e-10 and (k < epoch):
                k = k + 1
                if k == 1:
                    p = res.copy()
                else:
                    beta = rho/rho_h;p = res + beta*p
                w = np.zeros_like(res)
                for j in range(1,N - 1):
                    w[1:M - 1,j] = - r2*(p[1:M - 1,j - 1] + p[1:M - 1,j + 1]) - \
                    r1*(p[0:M - 2,j] + p[2:M,j]) + 2*(r1 + r2)*p[1:M - 1,j] 
                alpha = rho/(p*w).sum()
                u = u + alpha*p
                res = res - alpha*w;rho_h = rho;rho = (res*res).sum()
                #print(rho,np.linalg.norm(alpha*p))
            ela = time.time() - tic
            print('the end iteration:%d,the residual:%.3e,the time:%.2f'%(k,rho,ela),np.linalg.norm(alpha*p))
            return u,res
        def GS(self,u,rig,epoch):#rig也是一个矩阵,跟hx有关,hx控制rig的大小,扁平化处理以后就是方程的右端项
            #由于u已经初始化,把边界条件吸收进去,所以rig的边界为0
            tic = time.time()
            nx = [rig.shape[0],rig.shape[1]]
            hx = [(self.bound[0,1] - self.bound[0,0])/(nx[0] - 1),(self.bound[1,1] - self.bound[1,0])/(nx[1] - 1)]
            #print(rig.shape)
            u_new = u.copy()
            M = nx[0];N = nx[1]
            
            r1 = hx[1]/hx[0];r2 = hx[0]/hx[1]
            
            #GS
            
            for k in range(epoch): 
                for i in range(1,M - 1):
                    for j in range(1,N - 1):
                        u_new[i,j] = (hx[0]*hx[1]*rig[i,j] + r1*(u_new[i - 1,j] + u[i + 1,j]) + \
                        r2*(u_new[i,j - 1] + u[i,j + 1]))/(2*r1 + 2*r2)
                u = u_new.copy()
                #上面这部分迭代就是利用G-S迭代近似求解线性方程的过程,下面需要利用G-S迭代的结果求残差
            for j in range(1,N - 1):
                rig[1:M - 1,j] = hx[0]*hx[1]*rig[1:M - 1,j] - 2*(r1 + r2)*u_new[1:M - 1,j] + \
                r2*(u_new[1:M - 1,j - 1] + u_new[1:M - 1,j + 1]) + \
                r1*(u_new[0:M - 2,j] + u_new[2:M,j])
            ela = time.time() - tic
            print('the end iteration:%d,the time:%.2f'%(k,ela))
                
            return u_new,rig
            
                    
    bound = np.array([[0,4.0],[0,4.0]])
    
    prob = 1
    mg = MG(bound,prob)
    hx = [1/64,1/64]
    nx = [int((bound[0,1] - bound[0,0])/hx[0]) + 1,int((bound[1,1] - bound[1,0])/hx[1]) + 1]
    size = nx[0]*nx[1]
    X = np.zeros([size,2])
    m = 0
    for i in range(nx[0]):
        for j in range(nx[1]):
            X[m,0] = bound[0,0] + i*hx[0]
            X[m,1] = bound[1,0] + j*hx[1]
            m = m + 1
    u_acc = UU(X,[0,0],prob).reshape(-1,1)
    
    rig = FF(prob,X).reshape([nx[0],nx[1]])
    print(rig.shape)
    rig[0,:] = 0
    rig[-1,:] = 0
    rig[1:nx[0] - 1,0] = 0
    rig[1:nx[0] - 1,-1] = 0
    
    epoch = 2*size
    u_pred,res = mg.residual(mg.u_init(hx),rig,epoch)
    
    print(max(abs(u_acc - u_pred.reshape(-1,1))))
    
    
    

    在这里插入图片描述
    上面展示了共轭梯度法的结果,如果用高斯迭代,写完这篇博客还没跑出来。

    展开全文
  • 用超松弛迭代法求解接地金属槽内电位分布 一实验内容 试用超松弛迭代法求解接地金属槽内电位的分布 已知mm/4?a?4cm10h?a? ? 100 V 给定边值如图所示 给定初值)0?0?ji? 误差范围0?5?10?0? 计算迭代次数分布?ji? 0? ...
  • 非线性薄膜方程差分迭代求解,杨岳民,,本文对于类似于非线性薄膜一类的非线性偏微分方程,提出一种新的算法。该首先略去非线性部分的影响,或给予非线性部分某一初值
  • 用不同精度的差分格式将高维平稳 FPK方程离散化为线性代数方程组,然后用超松弛迭代法求解该线性代数方程组得到平稳 FPK方程的近似解。讨论了不同的差分格式、网格密度及超松弛因子对解精度及收敛速度的影响,并与其他...
  • 利用局部消元法建立了求解椭圆型方程的有限差分格式,并根据 Chebyshev多项式加速技术构造了一个混合半迭代法。该算法在第一层网格上仍使用经典的 Jacobi迭代法,在内层网格上使用多项式加速技术。数值实验表明,新算法...
  • 最近在学信号与系统,介绍了解差分方程的几种方法,迭代法、经典法和z变换法。迭代法用不上,但可以作为检验结果正确性的工具,还可用于求取经典法的初值。 想起来之前看过的一篇知乎文章,介绍利用卡西欧计算器实现...

    前言

    最近在学信号与系统,介绍了解差分方程的几种方法,迭代法、经典法和z变换法。迭代法用不上,但可以作为检验结果正确性的工具,还可用于求取经典法的初值。
    想起来之前看过的一篇知乎文章,介绍利用卡西欧计算器实现递归计算,所以粗略地研究了一下怎么用卡西欧解差分方程,其实利用的是最简单和基础的卡西欧编程技巧。这篇文章着实让我惊讶,因为利用卡西欧实现循环完全是利用卡西欧计算器的bug实现的,作者思路清奇。

    一个简单例子

    迭代求解: y(n)y(n1)y(n2)=0,y(0)=y(1)=1,n>=2y(n)-y(n-1)-y(n-2)=0,y(0)=y(1)=1,n>=2
    编程实现:先设初值,A=1,B=1,然后输入:
    ×AC:A+BA:CB\times A \rightarrow C: A+B \rightarrow A: C\rightarrow B
    很简单。这里记住循环n次,A=f(n+2)就行了。

    其他非齐次差分方程同理,在A+B这里添一项D,D迭代就行,一般都能解决。

    展开全文
  • 基于二维对流扩散方程的四阶紧致差分格式,将其转化为代数方程组,得到其三对角块形式的系数矩阵,利用稀疏矩阵存储技术和预条件迭代法进行求解,并与传统的中心差分格式所得数值解进行比较,充分说明了方法的高效性...
  • 用多重网格算法求解微分方程的matlab例子。程序采用采取四层网格,微分方程的离散选用有限差分法,每层网格上的计算采用逐次超松弛迭代法(SOR迭代);由细网格限制到粗网格,采用完全加权限制算子
  • 用有限差分法求解一维耦合偏微分方程组 通过极线等值线图显示磁场线的失真 对解的收敛速度和有限差分方法的稳定性进行数值分析 (最终会)演示如何使用 F2PY 极大地提高迭代代码的性能,同时处理 Python 中的所有更...
  • 描述:- 该程序的目标是求解区域 0<x<30, 0<y<30 中的稳态电压分布,假设正方形的一侧以 45*(x/xmax)* 的电压激发((1-x)/xmax)* 伏特 (xmax=30) 和所有其他边保持在 0 伏特。 该边界处的电压与其最大值...
  • 湖南理工学院学报.2014(03) 引言 线性常系数差分方程是描述线性时不变离散时间系统的数学模型 求解差分方程是分析离散时间系 统的重要内容.在信号与系统课程中介绍的求解方法主要有迭代法时域经典法双零法和变换域 ...
  • 这几天在家躲避疫情,闲来无事,写了这个多重网格法求解泊松方程的算法的代码。 多重网格法可能是目前为止解泊松方程最快的算法,...泊松方程化为差分方程后,每个格点都可以写成一个方程,因此得到一个方程组。使用...

    这几天在家躲避疫情,闲来无事,写了这个多重网格法求解泊松方程的算法的代码。
    多重网格法可能是目前为止解泊松方程最快的算法,n个格点需要n次计算就可以收敛,而快速傅里叶变换的收敛速度是n*logn, 共轭梯度法是n^2.。多重网格法可以方便的应对各种边界条件,这一点比傅里叶变换之类的谱方法要好得多。

    多重网格法可以这么理解。泊松方程化为差分方程后,每个格点都可以写成一个方程,因此得到一个方程组。使用迭代法解这个方程组的时候,粗网格可以快速消除低频误差,细网格消除高频误差,收敛速度远快于直接使用细网格做迭代。

    本程序写的时候只考虑了正方形网格的情形,但是可以很方便的改为长方形网格以及不等距长方形网格(目前不支持极坐标)。迭代方法采用了SOR松弛迭代法。例题选取的是我之前用共轭梯度法解的一个泊松方程。
    4.

    import numpy as np  
    import matplotlib.pyplot as plt
    import math
    from pylab import *  
    
    def fine2(grid0,hx0,hy0):   #化成细网格
        nx_grid2  =  int((hx0.size)*2 + 1) #x 方向粗化后网格数量
        ny_grid2  =  int((hy0.size)*2 + 1)        
        grid2  =  np.zeros((nx_grid2,ny_grid2),  dtype   =   float)
        hx2 = np.zeros((nx_grid2-1),  dtype   =   float)
        hy2 = np.zeros((ny_grid2-1),  dtype   =   float)
        x2 = np.zeros((nx_grid2),  dtype   =   float)
        y2 = np.zeros((ny_grid2),  dtype   =   float)
        x2[0] = xs
        y2[0] = ys
        for i in range(1, hx0.size+1):
            i2 = i*2
            for j in range(1, hy0.size+1):              #赋值中间的网格
                j2 = j*2
                grid2[i2,j2] = grid0[i,j]
                grid2[i2-1,j2] = grid0[i,j]*0.5+grid0[i-1,j]*0.5
                grid2[i2,j2-1] = grid0[i,j]*0.5+grid0[i,j-1]*0.5
                grid2[i2-1,j2-1] = grid0[i,j]*0.25+grid0[i,j-1]*0.25+grid0[i-1,j]*0.25+grid0[i-1,j-1]*0.25
        for i in range(0, hy0.size):#赋值边界
            grid2[0,2*i]  =  grid0[0,i]
            grid2[0,2*i+1]  =  0.5*grid0[0,i]+0.5*grid0[0,i+1] 
            grid2[-1,2*i]  =  grid0[-1,i]
            grid2[-1,2*i+1]  = 0.5* grid0[-1,i]+0.5*grid0[-1,i+1] 
        grid2[0,int(hy0.size*2)] =  grid0[0,hy0.size]  #边界最后一个点
        grid2[-1,int(hy0.size*2)] =  grid0[-1,hy0.size]  #边界最后一个点
    
        for j in range(0, hx0.size):#赋值边界
            grid2[2*j,0]  =  grid0[j,0]
            grid2[2*j+1,0]  =  0.5*grid0[j,0]+0.5*grid0[j+1,0] 
            grid2[2*j,-1]  =  grid0[j,-1]
            grid2[2*j+1,-1]  = 0.5* grid0[j,-1]+0.5*grid0[j+1,-1] 
        grid2[int(hx0.size*2),0] =  grid0[hx0.size,0]  #边界最后一个点
        grid2[int(hx0.size*2),-1] =  grid0[hx0.size,-1]  #边界最后一个点
    
        return grid2
        
    
    def hx_level(xsize):  #根据数组的大小判断是第几重网格,返回这一重的步长,x,y坐标
        level = int(math.log2((nx_grid-1)/xsize))#     计算粗化了几次
        #level_total  =  int(math.log2(phi.size-1))-5
        dstep  =  int(2**(level))
        hx2_grid  =  int((hx.size)/dstep)
        print(level,dstep,hx.size,hx2_grid)
        hy2_grid  =  int((hy.size)/dstep)
        hx2  =  np.zeros(hx2_grid,  dtype   =   float)
        hy2  =  np.zeros(hy2_grid,  dtype   =   float)
        for  i in range(0, hx2.size):                            
            j0  =  dstep*i
            h2_step  =  0.
            for j in range(dstep):    
                h2_step  =  h2_step + hx[j0 + j]
            hx2[i]  =  h2_step
        print(hx2)
    
        for  i in range(0, hy2.size):                            
            j0  =  dstep*i
            h2_step  =  0.
            for j in range(dstep):    
                h2_step  =  h2_step + hy[j0 + j]
            hy2[i]  =  h2_step
        x2 = np.zeros((hx2.size+1),  dtype   =   float)
        y2 = np.zeros((hy2.size+1),  dtype   =   float)
        x2[0] = xs
        y2[0] = ys
        for i in range(1, x2.size):# 赋值x2,默认不等距网格
            x2[i] = x2[i-1]+hx2[i-1]
        for j in range(1, y2.size):# 赋值y2,默认不等距网格
            y2[j] = y2[j-1]+hy2[j-1]
        return hx2,hy2,x2,y2
        
    
    def coarsen2(grid0,hx0,hy0):
        #dstep  =  int(2**(level-1))
        nx_grid2  =  int((hx0.size)/2 + 1) #x 方向粗化后网格数量
        ny_grid2  =  int((hy0.size)/2 + 1)
        grid2  =  np.zeros((nx_grid2,ny_grid2),  dtype   =   float)
        hx2 = np.zeros((nx_grid2-1),  dtype   =   float)
        hy2 = np.zeros((ny_grid2-1),  dtype   =   float)
        x2 = np.zeros((nx_grid2),  dtype   =   float)
        y2 = np.zeros((ny_grid2),  dtype   =   float)
        x2[0] = xs
        y2[0] = ys
        for i in range(1, nx_grid2-1):
            for j in range(1, ny_grid2-1):              #赋值中间的网格
               i2  =  2*i
               j2  =  2*j
               weight_total = 1./hx0[i2-1]+1./hx0[i2]+1./hy0[j2-1]+1./hy0[j2]
               weight1 = (1./hx0[i2-1])/weight_total     # 周边四个网格的权重
               weight2 = (1./hx0[i2])/weight_total
               weight3 = (1./hy0[j2-1])/weight_total
               weight4 = (1./hy0[j2])/weight_total
               grid2[i,j]  =  0.4*grid0[i2,j2] +  0.6*weight1*grid0[i2-1,j2] + 0.6*weight2*grid0[i2+1,j2]+0.6*weight3*grid0[i2,j2-1]+0.6*weight4*grid0[i2,j2+1]
            print(weight1,weight2,weight3,weight4)
        for i in range(0, ny_grid2):#赋值边界
            grid2[0,i]  =  grid0[0,2*i]
            grid2[-1,i]  =  grid0[-1,2*i]
        for j in range(0, nx_grid2):
            grid2[j,0]  =  grid0[2*j,0]
            grid2[j,-1]  =  grid0[2*j,-1]
    
        return grid2#,hx2,hy2,x2,y2      
    
    def Relax2(b2, phi0,  h,x1,x2):#注意,这个子程序还没来得及写成非均匀网格的,先用正方形网格凑合下
        omig = 1.85 #松弛银子
        leveli = int(math.log2((nx_grid-1)/h.size))
        print('level',leveli)
        ite  =  20*((leveli)*leveli)
        k = -4            #k 取负)值可保证至少迭代一定次数
        while(k <ite):
            print(k)
            for j in range(0,h.size+1):  #向上
                for i in range(1,h.size):       #沿着横坐标。      非均匀网格要改这里
                    
                    if j  =  =  0:                     #最低下一行
                        phi0[j,i] = np.copy(phi0[j+1,i])
                    elif j =  = h.size:
                        phi0[j,i] = np.copy(phi0[j-1,i])#最上一行
                    else:#五点差分法
                        phi0[j,i] = (1-omig)*np.copy(phi0[j,i])+omig*(np.copy(phi0[j,i-1])+np.copy(phi0[j,i+1])+np.copy(phi0[j-1,i])+np.copy(phi0[j+1,i])+h[i]*h[i]* b2[i,j])/4.  
            k = k+1
        return phi0
    
    def MG(b_mg,phi0, x1, x2, h):#多重网格法的子程序
        vh = Relax2(b_mg,phi0,  h, x1, x2)
        rh = residual(b_mg,vh, h, x1, x2)
    
        level_total  =  4#int(math.log2(phi.size-1))-4          #计算层数,V型,先粗化,再细化
        for i in range (1, level_total):   #coarse and iterate
            print(i)
            hx,hy,x11,x21 = hx_level(vh.shape[0]-1)
            hx2,hy2,x12,x22 = hx_level(int(hx.size/2))
            v2h0   =   coarsen2(vh,hx,hx)
            b2h  =  coarsen2(b_mg,hx,hx) 
            v2h = Relax2(b2h,v2h0,  hx2, x12, x22)
            vh   =   v2h
            b_mg  =  b2h
        for i in range (1, level_total):  # fine and itearate
            print(i)
            hx,hy,x11,x21 = hx_level(vh.shape[0]-1)
    
            hx2,hy2,x12,x22 = hx_level(int(hx.size*2))       
            b2h  =  fine2(b_mg,hx,hx) 
            v2h0   =   fine2(vh,hx,hx)
            v2h  =  Relax2(b2h,v2h0,  hx2, x12, x22)
            vh   =   v2h
            b_mg  =  b2h 
          
    
        return vh,x12,x22 
    
    
    
    
    
    def init(nx_grid,ny_grid,xs,xe,ys,ye,xleft_boundary_value,xright_boundary_value,dyleft_boundary_value,dyright_boundary_value):
        hx  =  (xe-xs)/(nx_grid-1)* np.ones(nx_grid-1,  dtype   =   float)   #x 方向步长,可以是非等距
        hy  =  (ye-ys)/(ny_grid-1)* np.ones(ny_grid-1,  dtype   =   float)   #y 方向步长
        x  =  np.zeros(nx_grid,  dtype   =   float)
        x[0]  =  xs
        x[-1]  =  xe
        for i in range(1, nx_grid-1):
            x[i]  =  x[i-1] + hx[i-1]
        y  =  np.zeros(ny_grid,  dtype   =   float)
        y[0]  =  ys
        y[-1]  =  ye
        for i in range(1, ny_grid-1):             # initialize x and y
            y[i]  =  y[i-1] + hy[i-1]
    
        phi   =   np.ones((x.size,y.size),  dtype   =   float)*0.1   # initialize results
        phi[:,0]  =  xleft_boundary_value
        phi[:,-1]  =  xright_boundary_value
    
        b  =  np.zeros((x.size,y.size),  dtype   =   float)#        initialize b
        for i in range(0, nx_grid):
            for j in range(0, ny_grid):
                b[i,j] = bxy(x[i],y[j])
        return x,y,hx,hy,phi,b
    
    def bxy(x0,y0):#方程右边的函数
        return -2.*math.pi*math.pi*math.sin(math.pi*x0)*math.cos(math.pi*y0)
    
    nx_grid  =  257 #at least 512 grids to reach enough depth
    xs  =  -1.
    xe  =  1.
    
    ny_grid  =  257  #at least 512 grids to reach enough depth
    ys  =  -1.
    ye  =  1.
    
    xleft_boundary_value = 0.     # u(x = -1) = u(x = 1) = 0
    xright_boundary_value = 0.
    dyleft_boundary_value = 0.    # du(y = -1) = u(y = 1) = 0
    dyright_boundary_value = 0.
    
    x,y,hx,hy,phi,b =  init(nx_grid,ny_grid,xs,xe,ys,ye,xleft_boundary_value,xright_boundary_value,dyleft_boundary_value,dyright_boundary_value)
    
    #grid2,hx2,hy2,x2,y2  = coarsen2(b,hx,hy)
    #grid3 = fine2(grid2,hx2,hy2)
    
    #print('hx',hx)
    #hxlevel2,hylevel2,x22,y22 = hx_level(hx2.size)
    
    phi2 = Relax2(b,phi,  hx,x,y)
    retult,x121,y121 = MG(b,phi, x, y, hx)
    
    plt.figure(1)
    print(retult.max(),retult.min())
    contourf(x121,  y121,  retult, 80,  cmap  =  'seismic')  
    plt.colorbar()
    '''
    plt.figure(2)
    contourf(y,  x,  phi2, 80,  cmap  =  'seismic')  
    
    plt.colorbar()
    '''
    
    plt.show()
    plt.close()
    

    结果如图
    在这里插入图片描述

    
    
    展开全文
  • 针对研究吊桥模型而建立的四阶微积分方程,提出利用有限差分法进行求解.采用Newton型迭代法处理非线性项,大大提高了收敛效率,并给出差分逼近的误差分析.数值算例说明了算法的可行性和有效性.
  • 在此代码中,泊松方程的边界条件是沿 4 个端壁的已知电位 100V 和 -100V。... 泊松方程使用有限差分法 (FDM) 迭代求解。 泊松方程的解被绘制为电势等值线。 电场使用梯度函数计算,也显示为颤动图。
  • 利用带填补数的不完全LU分解(ILUT(τ,s))作预处理器以及FGMRES(20)作迭代加速器,对非均匀网格上二维对流扩散方程的高精度紧致差分格式进行数值实验,并与均匀网格上的计算结果进行对比,数值结果显示出非均匀网格上...
  • 该代码采用有限差分格式来求解二维热方程。 位于任意值 1000 的计算域中心的加热块是初始条件。 底壁初始化为 100 个任意单位,是边界条件。 随着算法的推进,每 50 个时间步长使用一个电影函数来说明热扩散。 代码...
  • 利用传统方法很难在计算机上实现差分方程的解析解求解,提出了一种获得差分方程解析解的线性算法,该算法的基础是完全线形变化。其核心操作为降维处理,对高阶差分方程进行逐次降阶运算,直至获得其解析解表达式。...
  • 利用带填补数的不完全LU分解(ILUT(τ,s))作预处理器以及FGMRES(20)作迭代加速器,对非均匀网格上二维对流扩散方程的高精度紧致差分格式进行数值实验,并与均匀网格上的计算结果进行对比,数值结果显示出非均匀网格上...
  • 建立差分方程1.1 差分的定义1.2 差分方程2 差分方程的模拟框图2.1 基本部件单元2.2 由框图建立差分方程3 差分方程的经典解法3.1 递推迭代3.2 经典3.3 齐次解的常用函数形式3.4 特解的常用函数形式4 零输入响应的...
  • matlab求解偏微分方程,相比较把偏微分转成长分为方程再调用ode函数,利用离散差分法,使用迎风格式迭代求解数值解。
  • 递推方程求解方法

    万次阅读 多人点赞 2018-07-13 21:26:58
    1、迭代法不断用递推方程的右部替换左部,下面以汉诺塔为例进行求解。有时候直接迭代可能不太方便,可以使用换元迭代。下面以二归并排序迭代方程为例进行求解。2、消法 消法一般应用在递归方程右边不仅仅只...
  • 计算机控制系统

    2012-11-27 11:23:41
    计算机控制系统,实现实时监控、分时监控;增量式PID控制算法;递推迭代法求解差分方程
  • 递推方程求解的几种方法

    千次阅读 2020-10-08 17:17:50
    消法一般应用在递归方程右边不仅仅只依赖于当前项的前一项,而是前很多项,这种递归方程直接用迭代法很麻烦。属于高阶递归方程,因此要先把高阶递归方程进行消,再进行迭代。以快速排序的递归方程为例。 3、...
  • 使用迭代法求解 1:## 向前欧拉格式(注意,这种方法不稳定,不可取) ```python import numpy as np import matplotlib.pyplot as plt import torch import time import torch.nn as nn def UU(prob,X,t,ox,...
  • 用两同心圆将求解域划分为存在重叠的有限和无限两个区域,在有限和无限域上分别用有限元和差分线法求解Laplace方程边值问题。用差分线法推导出的关系式修正有限元方程,求解该方程组从而得到原问题的解。本算法将求解...
  • 目的 加速 SSOR迭代法的收敛性。方法 运用矩阵分裂理论及比较定理进行证明。结果 得到矩阵为严格对角占优 L-矩阵时,预条件后能够加速 ...结论 对于求解差分方法、有限元方法及科学计算中产生的线性方程组提供理论支持。
  • 利用降维法推导出非均匀网格上三维对流扩散方程的高精度紧致差分格式,对于离散得到的代数方程组采用BiCGStab(2)迭代法求解. 数值算例表明,在网格节点数相同的情况下,基于非均匀网格的计算格式较均匀网格格式...

空空如也

空空如也

1 2 3 4 5 ... 7
收藏数 128
精华内容 51
关键字:

迭代法求解差分方程