python中的高斯（-Legendre）求积答案

【问题标题】：Gauss(-Legendre) quadrature in pythonpython中的高斯（-Legendre）求积
【发布时间】：2015-01-22 19:39:45
【问题描述】：

我正在尝试使用高斯求积来近似函数的积分。（更多信息：http://austingwalters.com/gaussian-quadrature/）。第一个函数在区间 [-1,1] 上。第二个函数通过变量的变化推广到[a,b]。问题是我不断收到错误“'numpy.ndarray' object is not callable”。我假设（如果我错了，请纠正我）这意味着我试图将数组 w 和 x 调用为函数，但我不确定如何解决这个问题。

这是代码

from __future__ import division
from pylab import *
from scipy.special.orthogonal import p_roots

def gauss1(f,n):
    [x,w] = p_roots(n+1)
    f = (1-x**2)**0.5
    for i in range(n+1):
        G = sum(w[i]*f(x[i]))
    return G

def gauss(f,a,b,n):
    [x,w] = p_roots(n+1)
    f = (1-x**2)**0.5
    for i in range(n+1):
        G = 0.5*(b-a)*sum(w[i]*f(0.5*(b-a)*x[i]+ 0.5*(b+a)))
    return G

这些是相应的错误消息

gauss1(f,4)
Traceback (most recent call last):

  File "<ipython-input-82-43c8ecf7334a>", line 1, in <module>
    gauss1(f,4)

  File "C:/Users/Me/Desktop/hw8.py", line 16, in gauss1
    G = sum(w[i]*f(x[i]))

TypeError: 'numpy.ndarray' object is not callable

gauss(f,0,1,4)
Traceback (most recent call last):

  File "<ipython-input-83-5603d51e9206>", line 1, in <module>
    gauss(f,0,1,4)

  File "C:/Users/Me/Desktop/hw8.py", line 23, in gauss
    G = 0.5*(b-a)*sum(w[i]*f(0.5*(b-a)*x[i]+ 0.5*(b+a)))

TypeError: 'numpy.ndarray' object is not callable

【问题讨论】：

标签： python numpy gaussian numerical-integration

【解决方案1】：

在gauss1(f,n) 中，当f 是一个数组时，您将其视为一个函数，因为您正在重新分配它；

def gauss1(f,n):
  [x,w] = p_roots(n+1)
  f = (1-x**2)**0.5 # This line is your problem.
  for i in range(n+1):
    G = sum(w[i]*f(x[i]))
  return G

你在第二个函数中做了类似的事情。

【讨论】：

如果我将其更改为：def test(f,n): [x,w] = p_roots(n+1) for i in range(n+1): f[i] = (1-x[i]**2)**0.5 G = sum(w[i]*f([i])) return G 然后我得到错误：File "C:/Users/Me/Desktop/hw8.py", line 29, in test f[i] = (1-x[i]**2)**0.5: TypeError: 'builtin_function_or_method' object does not support item assignment:（抱歉缺少换行符......）
是的，那是因为您将一个函数传递给gauss1(f,n)，目的是创建一个近似它的高斯-勒让德正交。然后你应该在函数内部调用f(x) 来创建正交

【解决方案2】：

正如 Will 所说，您对数组和函数感到困惑。

需要单独定义要积分的函数，传入gauss。

例如

def my_f(x):
    return 2*x**2 - 3*x +15 

gauss(m_f,2,1,-1)

您也不需要循环，因为 numpy 数组是 vectorized 对象。

def gauss1(f,n):
    [x,w] = p_roots(n+1)
    G=sum(w*f(x))
    return G

def gauss(f,n,a,b):
    [x,w] = p_roots(n+1)
    G=0.5*(b-a)*sum(w*f(0.5*(b-a)*x+0.5*(b+a)))
    return G

【讨论】：

【解决方案3】：

quadpy，我的一个小项目，可能会有所帮助：

import numpy
import quadpy


def f(x):
    return numpy.exp(x)


scheme = quadpy.line_segment.gauss_legendre(10)
val = scheme.integrate(f, [0.0, 1.0])
print(val)

1.7182818284590464

一维还有许多其他的求积方案。

【讨论】：

【解决方案4】：

示例：使用 n = 2 的高斯积分求解 b = pi/2 和 a = 0 的积分 2+sinX

import numpy as np

E = np.array([-0.774597, 0.000000, 0.774597])
A = np.array([0.555556, 0.888889, 0.555556])

def gauss(f, a, b, E, A):
    x = np.zeros(3)
    for i in range(3):
        x[i] = (b+a)/2 + (b-a)/2 *E[i]

    return (b-a)/2 * (A[0]*f(x[0]) + A[1]*f(x[1]) + A[2]*f(x[2]))


f = lambda x: 2 + np.sin(x)
a = 0.0; b = np.pi/2

areaGau = gauss(f, a, b, E, A)
print("Gaussian integral: ", areaGau)

【讨论】：