多重网格法解泊松方程（两步法）

这里是一个一维泊松方程，采用两步法求解，教程见
http://bender.astro.sunysb.edu/classes/numerical_methods/lectures/elliptic-multigrid.pdf


import math
import numpy as np
import matplotlib.pyplot as pltdef coarsen(grid_fine,x_fine):# Lagrange type interpolation, two orderprint('xfine',x_fine)n_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):#   print(n_grid_fine ,i)i2 = np.count_nonzero(x_coarsen<x_fine[i])-1#find the index if x_fine[i] > x_coarsen[1]: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[i2])/(x_coarsen[i2]-x_coarsen[i2+1])l1  =  (x_fine[i]-x_coarsen[i2])/(x_coarsen[i2+1]-x_coarsen[i2])grid_fine[i]  =  l0*grid_coarsen[i2]+l1*grid_coarsen[i2+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  =  4**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 phi
def 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 r
def fbh(x):b = -16*np.sin(4*x) return b
def MG2(phi, x, b):print('x0',x)h0  =  x[1]-x[0]               vh  =  Relax2(b, phi, h0,2)rh  =  residual(b, vh, h0)# residualprint('x = ',x)i = 0while (np.max(rh)>1e-5): #print(phi)r2h,x2 =  coarsen(rh,x)   # residual to coarse gride2h0 = np.zeros(r2h.size,dtype = float)h2  =  x2[2]-x2[1]e2h  =  Relax2(r2h, e2h0,  h2,1)eh,x = fine(e2h,x2)b = fbh(x)v1h = vh+ehh = x[2]-x[1]phi = Relax2(b, v1h, h, 2)vh = v1hrh = residual(b, vh, h)i = i+1print(i)return phin_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#初始值，随便选，这里用了0.1
phi[0]  =  0.
phi[-1]  =  0.b  =  -16*np.sin(4*x)#泊松方程的右边，左边是一个phi的二阶导数
result2  =  np.sin(4*x)   #precise solution
result  =  MG2(phi, x, b)plt.figure(1)
plt.plot(x, result2, label  =  'original')
plt.scatter(x, result, label  =  'simulation')
plt.legend()
plt.show()
plt.close()

泊松方程的形式是phi’'(x)= -16np.sin(4x),精确解是phi(x)=sin(4*x),结果如图，可见精度还是比较高的

多重网格法解泊松方程（两步法）相关推荐

柱坐标下多重网格法解泊松方程-python
写了一个柱坐标下多重网格法解泊松方程的code,外边界采用的是第一类边界条件,直接用解析解赋值,内边界需要注意的是在x=0的转轴上面有奇点. 柱坐标下泊松方程形式为当x趋向于0时, 将方程离散化,其 ...
SAP-MM MM STO订单详解1（工厂间的转储一步法和两步法）
本文主要介绍一步法和两步法下的工厂间的转储,前提是这两个工厂隶属于同一个公司. 1.工厂间的转储:一步法案例:将物料30800807从工厂PG41通过一步法转移5PC至工厂PG42,操作之前需要判断 ...
量子计算机解泊松方程,试求泊松方程的解.ppt
试求泊松方程的解 2.5 具有非齐次边界条件的问题 (*) (14) (15) 利用公式其中系数满足那么其中系数计算可得 (94) 于是,问题(94)的解为因此,原问题(91)的解为求解 ...
Tsai两步法求手眼标定矩阵
Tsai两步法求手眼标定矩阵 Tsai方法介绍术语概念齐次变换矩阵和坐标系的定义旋转轴和旋转角度引理的证明和解释 AX=XB构造引理1 引理2 引理3 引理4 引理5 引理6 --得到公式( ...
模糊匹配读音_onenote搜索机制详解②：两种搜索模式，模糊与精确匹配
先从纯文本搜索讲起,这是最基本也是最重要的. 从这篇开始,以及接下来连续几篇文章,都会介绍搜索的基础功能.注意,这几篇文章中谈论的都是基本的.正常的搜索功能,暂时不考虑Bug等因素. 在很多软件(例如 ...
CV：利用cv2自带两步法haarcascade_frontalcatface.xml实现对猫脸检测
CV:利用cv2自带两步法haarcascade_frontalcatface.xml实现对猫脸检测目录输出结果实现代码输出结果实现代码 @author: niu ''' import cv ...
CV：利用cv2自带两步法haarcascade_frontalface_default.xml、_smile.xml实现对人脸、笑脸同时检测
CV:利用cv2自带两步法haarcascade_frontalface_default.xml._smile.xml实现对人脸.笑脸同时检测目录输出结果实现代码输出结果实现代码 #CV:利 ...
卡尔曼_卡尔曼估计两步法
在上一篇文章中手把手推导了一遍卡尔曼增益,不熟悉的小伙伴可以看养生的控制人:卡尔曼增益推导zhuanlan.zhihu.com 这里再回顾一下重点. 问题重述假设真实系统为其中 . 我们对系统 ...
用matlab给图片标记区域,MATLAB二值图像连通区域标记（两步法）
posted on2012-12-06 16:24Dsp Tian 两步法中第二步是比较麻烦的,其中用到了不相交集合的一些理论,尤其是不相交集合森林,这里的find_set函数就是参考<算法导论 ...

多重网格法解泊松方程（两步法）

多重网格法解泊松方程（两步法）相关推荐

最新文章

热门文章