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

2024-05-17 16:25:11 发布

您现在位置:Python中文网/ 问答频道 /正文

我正在尝试使用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.26816999999 四边形:-1598.494351085969

为什么结果不同?哪一个是正确的结果 这是函数的图表: enter image description here


Tags: 函数fromimportnumpyasnpfunctionscipy
1条回答
网友
1楼 · 发布于 2024-05-17 16:25:11

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

bad poly interpolation

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我们得到一个线性插值

如果你还想要别的东西,恐怕你得告诉我

相关问题 更多 >