python高斯求和函数_选择积分方法—高斯积分
无限区间 (1)梯形法则,(2)辛普森法则,(3)龙伯格积分法或(4)高斯积分法,有一些适用的指导原则。
通常,更高阶的方法对于平滑函数更好。 如果不是,那么使用更简单的方法会更好,因为数据的变化不会反映在采样点上。 梯形法则适用于在均匀间隔的采样点处积分来自实验的数据。 这对于表现不佳的函数是有好处的。 辛普森的规则依赖于被积函数的更高阶的近似,以便准确。 而高斯积分是非常准确的,如果你需要均匀间隔的采样点,它不是令人满意的。
高斯积分
######################################################################
#
# 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
# 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
# values of N up to 1000. It is compatible with version 2 and version
# 3 of Python.
#
# Written by Mark Newman , June 4, 2011
# You may use, share, or modify this file freely
#
##################%####################################################
from numpy import ones,copy,cos,tan,pi,linspace
def 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
计算实例
使用高斯积分的数值算法求解高斯积分
做积分变换
得到
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.
重写基本积分公式(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中的四点集进行比较。
映射积分点
GaussPoints.py
# GaussPoints.py: N point Gaussian quadrature pts & Wts generation
import numpy as np
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 = np.cos(np.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) )
python高斯求和函数_选择积分方法—高斯积分相关推荐
- python中bool函数用法_在python中bool函数的取值方法
bool是Boolean的缩写,只有真(True)和假(False)两种取值 bool函数只有一个参数,并根据这个参数的值返回真或者假. 1.当对数字使用bool函数时,0返回假(False),任何其 ...
- Python内置函数sorted()和列表方法sort()排序规则不得不说的事
Python内置函数sorted()和列表方法sort()可以使用key参数指定排序规则,并且都是稳定排序,也就是说,对于指定规则不能涵盖的元素,本来谁在前面,排好以后谁还是在前面. 直接用代码说话: ...
- python用circle函数画兔子的方法
python用circle函数画兔子的方法 circle函数说明 1.在circle函数中,参数radius取像素值和extent取角度的整数值可以取正负值. circle()函数以画笔当前方向(y' ...
- Python使用property函数定义属性访问方法如果不定义fget会怎么样?
我们知道Python使用property函数定义属性访问方法时的语法如下: 实例属性=property(fget=None, fset=None, fdel=None, doc=None) 而是要@p ...
- python高斯求和_利用Python进行数据分析(3)- 列表、元组、字典、集合
本文主要是对Python的数据结构进行了一个总结,常见的数据结构包含:列表list.元组tuple.字典dict和集合set. image 索引 左边0开始,右边-1开始 通过index()函数查看索 ...
- python帮助系统函数_【Python】【基础知识】【内置函数】【help的使用方法】
原英文帮助文档: help([object]) Invoke the built-in help system. (This function is intended for interactive ...
- python指数运算函数_分享Python中用于计算指数的exp()方法实例教程
exp()方法返回指数x: ex. 语法 以下是exp()方法的语法:import math math.exp( x ) 注意:此函数是无法直接访问的,所以我们需要导入math模块,然后需要用math ...
- python list find函数_对python中list的五种查找方法说明
Python中是有查找功能的,五种方式:in.not in.count.index,find 前两种方法是保留字,后两种方式是列表的方法. 下面以a_list = ['a','b','c','hell ...
- python列表去重函数_对python中两种列表元素去重函数性能的比较方法
测试函数: 第一种:list的set函数 第二种:{}.fromkeys().keys() 测试代码: #!/usr/bin/python #-*- coding:utf-8 -*- import t ...
- python中累加函数_对Python实现累加函数的方法详解
对Python实现累加函数的方法详解 发布时间:2020-10-26 00:02:44 来源:脚本之家 阅读:120 作者:岚漾忆雨 这个需求比较奇怪,要求实现Sum和MagaSum函数,实现以下功能 ...
最新文章
- java 动态获取类实例化_Java:使用反射动态实例化类
- python ndarray
- scala模式匹配match操作
- vbox下安装arch
- 微弱信号检测_世界上最轻薄的信号放大器:可精准监测生物信号!
- dell计算机运行慢怎么解决方法,戴尔笔记本电脑运行速度慢怎么办?
- 魔百和CM201-1-YS易视腾代工-免拆机-线刷固件包
- 超级表格企业版,最实用的三个功能
- 解决pip下载速度过慢及超时等其它的报错的方法适于多种操作系统(详细)
- Windows7双屏扩展及双屏桌面背景独立显示
- PS 模块BAPI新建修改项目、WBS、网络、作业 (一)
- google的秘密入口+搜索技巧
- [扩展阅读] EasyGUI 学习文档【超详细中文版】
- 自定义View基础之——canvas,paint的基本用法
- 关于圆周卷积和fft求卷积的一些看法
- 能用四川电信卡开通的虚拟服务器,双网通手机也能用电信卡了?VoLTE开放:发短信就能开通...
- solarwinds 快速上手指南
- python爬虫爬取百度搜索结果,Bob blog
- MATLAB 基础知识 数据类型 时间表 按行时间和变量类型选择时间表数据
- 兴牧生物科技获得波奇网千万级股权投资
热门文章
- ArcGIS单波段影像重分类与批处理
- 微型计算机原理考试试卷,微机原理试题集题库带答案
- 1047: 对数表 C语言
- linux scsi程序,Linux scsi设备读写流程
- 机器狗病毒样本 穿透冰和点还原卡
- 2021-02-18
- 如何从零开始系统运营微信公众号?
- amd显卡驱动目录linux,面向 Radeon、Radeon Pro、FirePro、APU、CPU、锐龙、台式机、笔记本的 AMD 驱动程序和支持...
- AHP计算权重.mat
- win10无法运行jre java_Windows10系统安装不了jre的解决方法