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

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

import math
import numpy as np  
import matplotlib.pyplot as plt

def coarsen(grid_fine,x_fine):# Lagrange type interpolation, two order
    print('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_fine
        l0  =  (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 interpolation

    grid_coarsen[0]  =  grid_fine[0]#left_boundary
    grid_coarsen[-1]  =  grid_fine[-1]#right_boundary
    return 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 interpolation
        else:
            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_boundary
    grid_fine[-1]  =  grid_coarsen[-1]#right_boundary

    return grid_fine, x_fine
    
def Relax2( b,  phi,  h,leveli):
    om  =  1.95
    ite  =  4**leveli    # the iteration increase with the level to get higher precision
    j   =   0
 
    while(j <ite): #控制迭代次数
        i   =   1                    # when the boundary is known,  set i   =   1
        while( 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 + 1
        j   =   j + 1
    return 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)# residual
    print('x = ',x)
    i = 0

    while (np.max(rh)>1e-5): 
        #print(phi)
        r2h,x2 =  coarsen(rh,x)   # residual to coarse grid
        e2h0 = 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+eh
        h = x[2]-x[1]
        phi = Relax2(b, v1h, h, 2)
        vh = v1h
        rh = residual(b, vh, h)
        i = i+1
    print(i)
    return phi


n_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),结果如图，可见精度还是比较高的