时间戳

首页 > 解决方案 > python scipy.integrate quad 和 trapz 给出了两个不同的结果

python - python scipy.integrate quad 和 trapz 给出了两个不同的结果

问题描述

我正在尝试使用 scipy 库集成一个函数：

import numpy as np
import pandas  as pd
from scipy.optimize import fsolve
from scipy.integrate import quad
from scipy.integrate import trapz
import matplotlib.pyplot as plt

x = np.array([-1,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1])*5.4
y = np.array([20.6398,  -45.2398,   -113.8779,  -52.7028,   618.7554,   -52.7028,   -113.8779,  -45.2398,   20.6398])
function = np.polyfit(x, y, 8)


def u(x):
  a = function[0]
  b = function[1]
  c = function[2]
  d = function[3]
  e = function[4]
  f = function[5]
  g = function[6]
  h = function[7]
  l = function[8]
  return  a*x**8 + b*x**7 + c*x**6 + d*x**5 + e*x**4 + f*x**3 + g*x**2 + h*x  + l

  
St =trapz(y,x)
print(St)

Squad, err = quad(u,-5.4,+5.4)
print(Squad)

结果是：trapz：291.2681699999999 quad：-1598.494351085969

为什么结果不同，哪个是正确的结果这是函数的图表： 在此处输入图像描述

标签： pythonnumpyscipy

解决方案

@Warren Weckesser 是正确的。您的多项式不能很好地插值。众所周知，高次多项式不适合这项任务。

t = np.linspace(x.min(), x.max(), 10**4)
plt.plot(t,u(t))
plt.plot(x,u(x))
plt.scatter(x,y)

坏多边形插值

trapz假设函数在您选择的点之间是线性的。所以如果你真的对这个多项式的积分感兴趣，你必须选择一个假设可以接受的步长。

t = np.linspace(x.min(), x.max(), 10**4)
trapz(u(t),t), quad(u,-5.4,+5.4)[0]

给你两次-1598.49... 我假设你对多项式不感兴趣，但我也不知道你想要什么。也许您想要图片中的线性插值？然后你可以这样做：

from scipy.interpolate import UnivariateSpline

spl = UnivariateSpline(x,y,k=1,s=0.1)
trapz(y,x), quad(spl,-5.4,+5.4)[0]

你得到两次 291.26...。请注意，如果您删除k=1并使用 s 参数UnivariateSpline进行插值会好得多。它也适用于多项式，但不是采用高次多项式来插入许多点，而是将多项式缝合在一起k。这就是k=1我们得到线性插值的原因。

如果你想要别的东西，恐怕你必须告诉我。

推荐阅读