howking的个人博客分享 http://blog.sciencenet.cn/u/howking

博文

计算物理——高斯积分代码

已有 2180 次阅读 2019-3-13 22:19 |个人分类:计算物理|系统分类:科研笔记




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.

image.png


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) )

image.png

计算结果


 Npts = 10 ,   Ans = 1.7182818284590295

 eps = 3e-14 , Error = -1.554312234475219e-14




https://blog.sciencenet.cn/blog-566204-1167377.html

上一篇:Monte Carlo 方法在统计物理中的应用 相关书籍和文献 yang元祐
下一篇:数据分析--拟合

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-12-3 20:23

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部