Gauss型求積公式
阿新 • • 發佈:2018-11-19
一、簡介
高斯求積公式是變步長數值積分的一種,基本形式是計算[-1,1]上的定積分。理論證明對於 n個節點的上述求積公式,最高有 2n - 1 次的代數精度,高斯公式就是使得上述公式具有 2n - 1次代數精度的積分公式。
二、實現
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 19 13:50:15 2016
Gauss型求積公式
@author: Administrator
"""
from numpy import *
from scipy import integrate
import matplotlib.pyplot as plt
#(-1,1)
def f1(x):
return x**2 / (1 - x**2)**(1/2)
#(0,pi/2]
def f2(x):
return sin(x) / x
def f2_symbol(x):
from sympy import sin
return sin(x) / x
def gauss_legendre(f,a,b,n):
from sympy import Symbol
from sympy import diff
from sympy import solve
from sympy import Eq
from math import factorial
x = Symbol('x')
l = diff((x**2 - 1)**(n+1),x,n+1)
l = l / (2**(n+1) * factorial(n+1))
root = solve(Eq(l,0))
A = []
result = 0
l_derivative = diff(l,x,1)
for i in range(0,n+1):
l_derivative_value = (l_derivative).evalf(subs={x:root[i]})
A.append(2 / ((1 - root[i]**2) * l_derivative_value**2))
result += (b -a) / 2 * A[i] * f((a+b)/2 + (b-a)/2*root[i])
return result
if __name__ == '__main__':
n = 4 # 1 2 4
print 'f1:'
print 'Truth-value:',integrate.quad(f1,-1,1)[0]
print 'Estimated-value:',float(gauss_legendre(f1,-1,1,n))
print
print 'f2:'
print 'Truth-value',integrate.quad(f2,0,pi/2)[0]
print 'Estimated-value:',float(gauss_legendre(f2_symbol,0,pi/2,n))
zero = 1e-10
x1 = linspace(-1,1,500)
x1[0] += zero
x1[-1] -= zero
y1 = f1(x1)
x2 = linspace(0,pi/2,500)
x2[0] += zero
y2 = f2(x2)
fig = plt.figure(figsize=(8,6))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax1.plot(x1,y1,color='b',linestyle='-',label='f1(x)')
ax1.legend(loc='upper left')
ax2.plot(x2,y2,color='r',linestyle='-',label='f2(x)')
ax2.legend(loc='lower left')
fig.show()
fig.savefig('a.png')