CODEWITHSUNDEEP
pythonnumpyscipy
NAS
NAS

Reputation: 218

python scipy.integrate quad and trapz gives two differents result

I am trying to integrate a function using the scipy library:

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)

the results are : trapz : 291.2681699999999 quad : -1598.494351085969

why the results are different and which is the correct result this is the graph of the function: enter image description here

Upvotes: 2

Views: 1714

Answers (1)

Lukas S
Lukas S

Reputation: 3603

@Warren Weckesser is correct. Your polynomial does not interpolate well. It is known that high degrees of polynomials are not suited for this task.

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 assumes that the function is is ruffly linear between the points you choose. So if you're actually interested in the integral of this polynomial you have to pick a step size for which the assumption is ok.

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

gives you twice -1598.49... I assume you're not interested in the polynomial but I also don't know what you want. Maybe you wanted the linear interpolation as in your picture? Then you could do:

from scipy.interpolate import UnivariateSpline

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

And you get twice 291.26.... Notice that if you remove the k=1 and play around with the s parameter UnivariateSpline is much much better for interpolation. It also work with polynomials but instead of taking a high degree polynomial to interpolate many points it stitches together polynomials of degree k. That's why for k=1 we get a linear interpolation.

And if you want something else I am afraid you'll have to tell me.

Upvotes: 1

Related Questions