||
import math from numpy import * def GaussPoints(Npts, a, b, x, w, eps): m = 0; i = 0; j = 0; t = 0.; t1 = 0.; pp = 0. p1 = 0.; p2 = 0.; p3 = 0. m = int((Npts+1)/2) for i in range(1, m+1): t = math.cos(math.pi*(float(i)-0.25)/(float(Npts)+0.5)) t1 = 1 while((abs(t-t1)) >= eps): p1 = 1. ; p2 = 0. for j in range(1, Npts + 1): p3 = p2; p2 = p1 p1 = ((2.*float(j)-1)*t*p2 - (float(j)-1.)*p3)/(float(j)) pp = Npts*(t*p1 - p2)/(t*t - 1.) t1 = t t = t1 - p1/pp x[i-1] = -t x[Npts-i] = t w[i-1] = 2./( (1.-t*t)*pp*pp) w[Npts-i] = w[i-1] for j in range(0, Npts): # Scale [-1,+1] to [a,b] x[j] = x[j]*(b-a)/2. + (b+a)/2. w[j] = w[j]*(b-a)/2.
from numpy import *; from GaussPoints import GaussPoints Npts = 10; Ans = 0; a = 0.; b = 1.; eps = 3.E-14 w = zeros(2001, float); x = zeros(2001, float) # Arrays def f(x): return exp(x) # Integrand GaussPoints(Npts, a, b, x, w, eps) # eps: precison of pts for i in range(0,Npts): Ans += f(x[i])*w[i] # Sum integrands print ('\n Npts =', Npts, ', Ans =', Ans) print (' eps =',eps, ', Error =', Ans-(exp(1)-1) )
计算结果
Npts = 10 , Ans = 1.7182818284590295
eps = 3e-14 , Error = -1.554312234475219e-14
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-12-26 06:03
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社