多重网格法（multigrid）求解1d泊松方程--python

多重网格法（multigrid）求解1d泊松方程

多重网格法可能是目前已知最快的求解泊松方程的算法，特别是在边界复杂的情况下，实用性远高于快速傅里叶算法。我写里写了一个一维求解泊松方程的程序,供大家参考。里面使用了三点式拉尔朗日插值进行粗细网格之间的转化，松弛迭代法解泊松方程。泊松方程的系数矩阵是对角线为-2，对角线两侧为1，边界为u0（已知），未知数从u1开始。多重网格采用了V型结构，共三层，外面套了一个W型循环。

教程参考
http://bender.astro.sunysb.edu/classes/numerical_methods/lectures/elliptic-multigrid.pdf
http://staff.ustc.edu.cn/~humaobin/course/cht/ppt/8.7.pdf
直接上代码


import math
import numpy as np
import matplotlib.pyplot as pltdef coarsen(grid_fine,x_fine):# Lagrange type interpolation, two ordern_grid_coarsen  =  int(grid_fine.size/2)grid_coarsen  =  np.zeros(n_grid_coarsen,  dtype   =   float)x_coarsen = np.linspace(x_fine[0],x_fine[-1],n_grid_coarsen)for i in range(1, grid_coarsen.size):# print(grid_coarsen.size, i)i2  = np.count_nonzero(x_fine<x_coarsen[i])-1# find the index in x_finel0  =  (x_coarsen[i]-x_fine[i2])*(x_coarsen[i]-x_fine[i2+1])/(x_fine[i2-1]-x_fine[i2])/(x_fine[i2-1]-x_fine[i2+1])l1  =  (x_coarsen[i]-x_fine[i2-1])*(x_coarsen[i]-x_fine[i2+1])/(x_fine[i2]-x_fine[i2-1])/(x_fine[i2]-x_fine[i2+1])l2  =  (x_coarsen[i]-x_fine[i2-1])*(x_coarsen[i]-x_fine[i2])/(x_fine[i2+1]-x_fine[i2-1])/(x_fine[i2+1]-x_fine[i2])grid_coarsen[i]  =  l0*grid_fine[i2-1] + l1*grid_fine[i2] +l2*grid_fine[i2 + 1]#Lagrange type interpolationgrid_coarsen[0]  =  grid_fine[0]#left_boundarygrid_coarsen[-1]  =  grid_fine[-1]#right_boundaryreturn grid_coarsen, x_coarsen
def fine(grid_coarsen,x_coarsen): #n_grid_fine  =  int(grid_coarsen.size*2)grid_fine  =  np.zeros(n_grid_fine,  dtype   =   float)x_fine  =  np.linspace(x_coarsen[0], x_coarsen[-1],n_grid_fine)for i in range(1, grid_fine.size):if x_fine[i] > x_coarsen[1]:i2 = np.count_nonzero(x_coarsen<x_fine[i])-1#find the index l0  =  (x_fine[i]-x_coarsen[i2])*(x_fine[i]-x_coarsen[i2+1])/(x_coarsen[i2-1]-x_coarsen[i2])/(x_coarsen[i2-1]-x_coarsen[i2+1])l1  =  (x_fine[i]-x_coarsen[i2-1])*(x_fine[i]-x_coarsen[i2+1])/(x_coarsen[i2]-x_coarsen[i2-1])/(x_coarsen[i2]-x_coarsen[i2+1])l2  =  (x_fine[i]-x_coarsen[i2-1])*(x_fine[i]-x_coarsen[i2])/(x_coarsen[i2+1]-x_coarsen[i2-1])/(x_coarsen[i2+1]-x_coarsen[i2])grid_fine[i]  =  l0*grid_coarsen[i2-1] + l1*grid_coarsen[i2] +l2*grid_coarsen[i2 + 1]#Lagrange type interpolationelse:l0  =  (x_fine[i]-x_coarsen[1])/(x_coarsen[0]-x_coarsen[1])l1  =  (x_fine[i]-x_coarsen[0])/(x_coarsen[1]-x_coarsen[0])grid_fine[i]  =  l0*grid_coarsen[0]+l1*grid_coarsen[1]grid_fine[0]  =  grid_coarsen[0]#left_boundarygrid_fine[-1]  =  grid_coarsen[-1]#right_boundaryreturn grid_fine, x_finedef Relax2( b,  phi,  h,leveli):om  =  1.95ite  =  3**leveli    # the iteration increase with the level to get higher precisionj   =   0while(j <ite): #控制迭代次数i   =   1                    # when the boundary is known,  set i   =   1while( i < phi.size-1):  phi_i  =  np.copy(phi[i])phi[i]  =  (1.-om)*phi_i+om*0.5*(phi[i + 1] + phi[i-1]-(h*h)*b[i])i   =   i + 1j   =   j + 1return phidef residual(f, u, h):r  =  np.zeros(f.size,  dtype   =   float)for i in range(1, u.size-1):r[i]  =  f[i]*(h*h)-(u[i + 1]-2.*u[i] + u[i-1])return rdef fbh(x):b = -64*np.sin(8*x) return b
def MG2(phi, x, b):h0  =  x[1]-x[0]               vh0  =  Relax2(b, phi, h0,1)rh  =  residual(b, vh0, h0)# residualeh0 = np.zeros(rh.size,dtype = float)eh =  Relax2(rh, eh0,  h0,1)vh = vh0+ehrh  =  residual(b, vh, h0)t = 0while ((np.max(rh) > 1e-2) and (t < 10)): for i in range(1,5): k = ir2h,x2 =  coarsen(rh,x)   # residual to fine gride2h0 = np.zeros(r2h.size,dtype = float)h2  =  x2[2]-x2[1]e2h  =  Relax2(r2h, e2h0,  h2,k)print(vh.size,x.size)v2h,x2 = coarsen(vh,x)phi = v2h+e2hb = fbh(x2)vh =  Relax2(b, phi, h2,k)rh = residual(b, vh, h2)x = x2for i in range(1,5): k = 6-ir2h,x2 =  fine(rh,x)   # residual to fine gride2h0 = np.zeros(r2h.size,dtype = float)# initiate the errh2  =  x2[2]-x2[1]e2h  =  Relax2(r2h, e2h0,  h2,k)     #calculate errprint(vh.size,x.size)v2h,x2 = fine(vh,x)                  #to fine gridphi = v2h+e2h                        # update phib = fbh(x2)                          # calculate b on this levelvh =  Relax2(b, phi, h2,k)           # calculate  vh again with updated phirh = residual(b, vh, h2)             #calculate residualx = x2t = t+1print(t)return vh,x2n_grid  =  128#should be 2**n
xs  =  0.0
xe  =  math.pi
x = np.linspace(xs,xe,n_grid)
h = x[1]-x[0]
phi   =   np.ones(x.size,  dtype   =   float)*0.1
phi[0]  =  0.
phi[-1]  =  0.b  =  -64*np.sin(8*x)#泊松方程，b是二阶导数的值。泊松方程形式phi''=b#precise solution
result,x2  =  MG2(phi, x, b)
result2  =  np.sin(8*x2) plt.figure(1)
plt.plot(x2, result2, label  =  'original')
plt.scatter(x2, result, label  =  'simulation')
plt.legend()
plt.show()
plt.close()

结果如图

多重网格法（multigrid）求解1d泊松方程--python相关推荐

多重网格法-松弛迭代法-二维泊松方程-python实现
这几天在家躲避疫情,闲来无事,写了这个多重网格法求解泊松方程的算法的代码. 多重网格法可能是目前为止解泊松方程最快的算法,n个格点需要n次计算就可以收敛,而快速傅里叶变换的收敛速度是n*logn, 共 ...
柱坐标下多重网格法解泊松方程-python
写了一个柱坐标下多重网格法解泊松方程的code,外边界采用的是第一类边界条件,直接用解析解赋值,内边界需要注意的是在x=0的转轴上面有奇点. 柱坐标下泊松方程形式为当x趋向于0时, 将方程离散化,其 ...
em算法python代码_EM 算法求解高斯混合模型python实现
注:本文是对<统计学习方法>EM算法的一个简单总结. 1. 什么是EM算法? 引用书上的话: 概率模型有时既含有观测变量,又含有隐变量或者潜在变量.如果概率模型的变量都是观测变量,可以直接 ...
python查询斐波那契数列通项公式_斐波那契数列求解总结（Python版）
最近在查阅斐波那契数列时,看到下面的文章,总结得非常好,于是自己上手使用 Python 练习并实现多种求解方法守望:面试官问你斐波那契数列的时候不要高兴得太早zhuanlan.zhihu.com ...
计算机软件求解线性规划模型--Python
线性规划简介线性规划(Linear programming,简称LP),是运筹学中研究较早.发展较快.应用广泛.方法较成熟的一个重要分支,它是辅助人们进行科学管理的一种数学方法.研究线性约束条件下线 ...
规划求解 python_使用Python/PuLp解决线性规划问题
Python中有许多第三方的工具可以解决这类问题,本教程介绍pulp工具包的使用. 一 · 教程讲解视频二 · 常用的线性规划求解软件 1.Excel 2.Lingo 3.Matlab 三 · Py ...
机器学习之利用SMO算法求解支持向量机—基于python
大家好,我是带我去滑雪! 本期将讨论支持向量机的实现问题,我们知道支持向量机的学习问题可以化为求解凸二次规划问题.这样的凸二次规划问题具有全局最优解,并且有许多最优化算法可以用于这一问题的求解.但是当 ...
期末微积分考试试题求解：利用python求解
在今年期末微积分考试试题:看看你能够在两个小时内做对几道题? 搜集到了一份期末微积分考试试题.为了对其内容进行进一步分析,对其内容进行整理如下. §01 填空题每个空3分,共10题 1. 求解常微 ...
牛顿法求解方程（python和C++）
牛顿法 GitHub: https://github.com/Sean16SYSU/Algorithms4N 目标函数为: f(x)=0f(x) = 0f(x)=0 算法很简单如下: xn+1=xn− ...

多重网格法（multigrid）求解1d泊松方程--python

多重网格法（multigrid）求解1d泊松方程

多重网格法（multigrid）求解1d泊松方程--python相关推荐

最新文章

热门文章