二重积分的复化Simpson方法

一般数值计算教材都只讲一重积分，很少讲二重积分和三重积分，这里我写了一个用Simpson方法求解二重积分的程序，供大家参考。程序的要点在于构建一个系数矩阵，具体教程可以参考
https://max.book118.com/html/2018/0419/162077258.shtm（张正印，二重积分的Simpson公式及其误差估计）

import math
import numpy as np  
import matplotlib.pyplot as plt
import sys
# this code calculate the 2d intergration with simpson method, the coffecient matrix is fixed and being created first, then the results can easily be calculated. . 

def grid_line1(n):# n is the size of x, should be odd
    grid=np.zeros(n,dtype=float)
    for i in range (0,n):
        if i%2==0:
            grid[i]=2
        else:
            grid[i]=4
    grid[0]=grid[-1]=1
    return grid
def grid_line2(n):# n is the size of x, should be odd
    grid=np.zeros(n,dtype=float)
    for i in range (0,n):
        if i%2==0:
            grid[i]=8
        else:
            grid[i]=16
    grid[0]=grid[-1]=4
    return grid
def grid_line3(n):# n is the size of x, should be odd
    grid=np.zeros(n,dtype=float)
    for i in range (0,n):
        if i%2==0:
            grid[i]=4
        else:
            grid[i]=8
    grid[0]=grid[-1]=2
    return grid
def simpson2d_grid(x,y):
    n=x.size
    m=y.size
    if n%2==0 or m%2==0:
        print("error, the size of x should be odd")
        sys.exit()
    else:
        grid1=grid_line1(n)
        grid2=grid_line2(n)
        grid3=grid_line3(n)
        grid0 = np.zeros((m,n),dtype=float)

        for i in range (0,m):
            if  i%2==0:
                grid0[i,:]=np.copy(grid3)

            else:
                grid0[i,:]=np.copy(grid2)
        grid0[0,:]=np.copy(grid1)
        grid0[-1,:]=np.copy(grid1)

    return grid0
def b(x0,y0):# the function is intergrate2d(ln(x+2*y)dxdy),x(1.4-2.0),y(1.0-1.5)

    return math.log(x0+2*y0)

xs = 1.4
xe = 2.0
ys = 1.0
ye = 1.5
n_x_grid=9
n_y_grid=9

x=np.linspace(xs,xe,n_x_grid)
y=np.linspace(ys,ye,n_y_grid)

coef_matrix=simpson2d_grid(x,y)
print(coef_matrix)
h=x[1]-x[0]
k=y[1]-y[0]
result = 0

for j in range(0,y.size):
    for i in range(x.size):
        result=result+h*k/9*coef_matrix[j,i]*b(x[i],y[j])
print(result)

系数矩阵和结果如下