>>>>>>>>>>>>>> ## Functions to calculate integration points and weights for Gaussian
# quadrature
#
# x,w = gaussxw(N) returns integration points x and integration
# weights w such that sum_i w[i]*f(x[i]) is the Nth-order
秀品网# Gaussian approximation to the integral int_{-1}^1 f(x) dx
# x,w = gaussxwab(N,a,b) returns integration points and weights
# mapped to the interval [a,b], so that sum_i w[i]*f(x[i])
# is the Nth-order Gaussian approximation to the integral
# int_a^b f(x) dx
#
# This code finds the zeros of the nth Legendre polynomial using
# Newton's method, starting from the approximation given in Abramowitz
动易2006# and Stegun 22.16.6. The Legendre polynomial itself is evaluated
# using the recurrence relation given in Abramowitz and Stegun
# 22.7.10. The function has been checked against other sources for
machine civilization# values of N up to 1000. It is compatible with version 2 and version
# 3 of Python.
#
# Written by Mark Newman <mejn@umich.edu>, June 4, 2011
# You may use, share, or modify this file freely
#
>>>###%>>>>>>>>>>## from numpy import ones,copy,cos,tan,pi,linspacedef gaussxw(N):
# Initial approximation to roots of the Legendre polynomial
a = linspace(3,4*N-1,N)/(4*N+2)
x = cos(pi*a+1/(8*N*N*tan(a)))
# Find roots using Newton's method
epsilon = 1e-15
delta = 1.0
while delta>epsilon:
p0 = ones(N,float)
p1 = copy(x)
追风伞for k in range(1,N):
p0,p1 = p1,((2*k+1)*x*p1-k*p0)/(k+1)
dp = (N+1)*(p0-x*p1)/(1-x*x)
dx = p1/dp
x -= dx
delta = max(abs(dx))
# Calculate the weights
w = 2*(N+1)*(N+1)/(N*N*(1-x*x)*dp*dp)
return x,w
def gaussxwab(N,a,b):
x,w = gaussxw(N)
return 0.5*(b-a)*x+0.5*(b+a),0.5*(b-a)*w
计算实例
使⽤⾼斯积分的数值算法求解⾼斯积分
做积分变换
best jeanist得到
from gaussxw import gaussxwab
from math import exp
def f(z):
return exp(-z**2/(1-z)**2)/(1-z)**2
N = 50
a = 0.0
b = 1.0
x,w = gaussxwab(N,a,b)
s = 0.0
for k in range(N):
s += w[k]*f(x[k])
print(s)
the value of the actual integral is √π/2, and Gaussian quadrature is accurate to machine precision.
Gaussian quadrature
重写基本积分公式（6.3）通常很有⽤，这样我们就可以将加权函数W（x）与被积函数分开：
⾼斯求积⽅法, N个点和权重选择的近似误差消失如果g (x)是⼀个(2 N-1)度多项式。为了得到这个不可思议的优化，点 最终在[a, b]上有⼀个特定的分布。⼀般情况下，如果g(x)是光滑的，或者可以通过提出⼀些W(x)来使其光滑(表6.2.4)，那么对于相同数量的点，⾼斯算法将⽐梯形规则和⾟普森规则产⽣更⾼的精度。有时被积函数可能不是光滑的，因为它在不同的区域有不同的⾏为。在这些情况下，分别对每个区域进⾏积分，然后将答案相加是有意义的。事实上，⼀些“智能”积分⼦例程决定使⽤多少个区间以及在每个区间中
使⽤什么规则。
表6.2.4所⽰规则均为⾼斯分布，⼀般形式为(6.32)。我们可以看到，在⼀种情况下权重函数是指数函数，在另⼀种情况下是⾼斯函数，在⼏种情况下是可积分奇点。与等间距规则不同，在区间的极值处从来不存在积分点，⽽点的值和权值随着点的个数N的变化⽽变化。虽然我们将离开⾼斯点的推导和权重数值⽅法的引⽤，我们注意，对于普通⾼斯(Gauss-Legendre)积分，然后点易变成零的勒让德多项式，权重0相关的微分,
。在数学函数库中，⽣成这些点和权重的⼦例程是标准的，可以在[A&S 72]之类的表中到，或者可以计算。我们提供的⾼斯⼦例程还将点缩放到指定区域。为了检查您的观点是否正确，您可能需要将它们与表6.1中的四点集进⾏⽐较。
映射积分点
中国麻风皮肤病杂志
Legendre polynomials勒让德多项式
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议