谈谈高维积分的数值方法 用Python给我写几个例子

高维积分数值方法是处理高维度函数的一种方法。通常，计算高维积分的解析解是非常困难的，因此我们需要使用数值方法来求解。以下是一些在Python中实现高维积分数值方法的例子。

1. 蒙特卡罗方法

蒙特卡罗方法是一种基于随机采样的数值方法。对于一个函数 f(x1, x2, ..., xn) 在超立方体 [a1, b1] x [a2, b2] x ... x [an, bn] 上的 n 维积分可以通过采样 N 个点的方法计算。其中每个采样点 (x1i, x2i, ..., xni) 都是独立且均匀分布在超立方体上，然后使用如下公式求积分值 I：

I = V * (1/N) * sum( f(xi1, xi2, ..., xin) )

这里，V 是超立方体的体积。

以下是 Python 代码实现：

import numpy as np

def monte_carlo_integration(f, limits, N):
V = np.prod([b-a for a, b in limits]) # 计算超立方体的体积
I = 0
for _ in range(N):
x = [np.random.uniform(a, b) for a, b in limits]
I += f(*x)
I *= V / N
return I

示例：

# 计算高维正态分布的积分
from scipy.stats import multivariate_normal

p = multivariate_normal(mean=[0,0,0], cov=np.identity(3)) # 三维正态分布

I = monte_carlo_integration(p.pdf, [(-1, 1), (-1, 1), (-1, 1)], int(1e6))
print(I)

2. 复合梯形法

复合梯形法是一种基于插值的数值方法。它将积分区间分为若干个子区间，并且在每个子区间上使用梯形公式计算积分。然后将这些积分值相加，得到最终的积分值。

以下是 Python 代码实现：

def composite_trapezoidal_rule(f, limits, N):
a, b = limits
h = (b-a) / N
x = np.linspace(a, b, N+1)
y = f(x)
I = h * (y[0] + y[-1] + 2 * np.sum(y[1:-1]))
return I

示例：

# 计算二维函数的积分
def f(x, y):
return np.sin(x) * np.cos(y)

I = composite_trapezoidal_rule(lambda x: composite_trapezoidal_rule(lambda y: f(x, y), [0, 1], 100), [0, 1], 100)
print(I)

3. 复合 Simpson 法

复合 Simpson 法是一种基于插值的数值方法。它将积分区间分为若干个子区间，并且在每个子区间上使用 Simpson 公式计算积分。然后将这些积分值相加，得到最终的积分值。

以下是 Python 代码实现：

def composite_simpsons_rule(f, limits, N):
a, b = limits
h = (b-a) / N
x = np.linspace(a, b, N+1)
y = f(x)
I = h/3 * (y[0] + y[-1] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-1:2]))
return I

示例：

# 计算三维函数的积分
def f(x, y, z):
return x**2 * y**3 * z**4

I = composite_simpsons_rule(lambda x: composite_simpsons_rule(lambda y: composite_simpsons_rule(lambda z: f(x, y, z), [0, 1], 10), [0, 1], 10), [0, 1], 10)
print(I)