对于复杂的积分怎么计算呢?这就提出了数值积分,他就是通过将连续区间进行离散,对离散量求和便是积分的近似值。
1D积分 ∫fdξ 在-1到1上求积分,数值积分就是,将-1到1离散为n个点,ΣAk*f(ξk) 其中ξk为-1,-1+2/n,-1+4/n,....,1。Ak为对应的ξk在每个点处的权系数。
=========================================
高斯积分:数值积分的一种。规定了Ak与ξk怎么求。
下面展开叙述,如何求Ak与ξk:
1. 1D积分 ∫fdx 在-1到1上求积分 我下面的积分区间都是-1到1,没写懒得打了。
1.1 1点高斯积分
∫fdξ = A1*f(ξ1) 式1
令 f=a0 则式1化为 2*a0 = A1*a0 式子1.1
令 f=a1*ξ 则式1化为 0 = A1*a1*ξ1 式子1.2
通过式1.1与1.2可求出 A1=2 ξ1=0
所以1点高斯积分为 2*f(0)
1.2 2点高斯积分
∫fdξ = A1*f(ξ1) + A2*f(ξ2) 式1
令 f=a0 则式1化为 2*a0 = A1*a0 + A2*a0 式子1.1
令 f=a1*ξ 则式1化为 0 = A1*a1*ξ1 + A2*a1*ξ1 式子1.2
令 f=a2*ξ**2 则式1化为 0 = A1*a2*ξ2**2 + A2*a2*ξ2**2 式子1.3
令 f=a3*ξ**3 则式1化为 0 = A1*a3*ξ3**3 + A2*a3*ξ3**3 式子1.4
通过式1.1与1.2 1.3 1.4, 四个式子四个未知量,可求出 A1=1 A2=1 ξ1=1/sqrt(3) ξ2=1/sqrt(3)
所以2点高斯积分为 1*f(-1/sqrt(3)) + 1*f(1/sqrt(3))
2. 2D积分 ∫∫fdxdy x为-1到1 y为-1到1
下面就举个例子如何利用2点高斯积分求解2D积分。
f = sin(exp(x*y))
先以y为常数,则∫fdx = sin(exp(-1/sqrt(3)*y)) + sin(exp(-1/sqrt(3)*y))
接着再对sin(exp(-1/sqrt(3)*y)) 求积分 ∫sin(exp(-1/sqrt(3)*y)) dy = sin(exp(-1/3)) + sin(exp(1/3))
接着对sin(exp(1/sqrt(3)*y))求积分 ∫sin(exp(1/sqrt(3)*y)) dy = sin(exp(1/3)) + sin(exp(-1/3))
则∫∫fdxdy (x为-1到1 y为-1到1) 该2D的两点高斯积分便是 2*sin(exp(-1/3)) + 2*sin(exp(1/3))
即2D与3D的高斯积分都可以转化为1D的高斯积分计算,只需要求出1D积分中的积分点与权系数即可。一般也不用求,查表即可。
=============================================
补充一下,自己用python算的结果。
不知道为什么sympy不出结果,只是一个表达式,所以我用scipy计算二重积分。
代码
import sympy as sp
x, y = sp.symbols('x, y')
# sympy计算结果
val1 = sp.integrate(sp.sin(sp.exp(x*y)),(x,-1,1), (y,-1,1))
# =======================================
from scipy import integrate
import numpy as np
val2,err2 = integrate.dblquad(lambda y,x:np.sin(np.exp(x*y)),-1,1,-1,1) # 参数依次为表达式 x下限 x上限 y下限 y上限
# scipy计算结果
print ('二重积分结果:',val2)
# ===================================
# 高斯积分
val1 = 2*(np.sin(np.exp(1/3)) + np.sin(np.exp(-1/3)))
print('高斯积分结果', val1)