python - 如何使用 sympy 找出定积分近似值的误差？

Question

我们的任务是使用 sympy 计算函数的精确定积分，然后将其与我们编写的另一个 python 函数获得的定积分近似值进行比较。使用简单的多项式函数，我的代码可以正常工作，但是使用复杂的正弦函数，它会不断中断或返回 nan。

from numpy import *
def simpson(f,a,b,n):
  if n<=0 or n%2!=0:
    print('Error: the number of subintervals must be a positive even number')
    return float('NaN')
  h = float(b - a) / float(n)
  x = arange(a,b+h,h)
  fx = f(x)
  fx[1:n:2] *= 4.0
  fx[2:n:2] *= 2.0
  return (h/3.)*sum(fx)

这是在一个文件 (simpandtrap) 中，使用具有 n 个子区间的辛普森法则近似给出 f 从 a 到 b 的定积分的近似

from pylab import *

def s(x):
  return x*sin(3./(x+(x==0)))    

这是给我带来麻烦的函数，在一个名为 assignment8functions 的文件中

import assignment8functions as a
import SimpAndTrap as st
import sympy as sp
x = sp.symbols('x')
Exact_int_q = sp.integrate(a.q(x),(x,0,2)).evalf(25)
Exact_int_s = sp.integrate(x*sp.sin(3./(x)),(x,0,2)).evalf(25)

q(x) 是我们应该使用的另一个函数，它一切正常 - 它只是一个多项式。当我尝试以相同的方式进行集成时，它会中断，所以我不得不将 s(x) 的函数直接放入调用中，而不是从另一个文件导入函数

n = a.array([10,100,1000,10000,10000,1000000])
s_error_simp_array = a.zeros(6)
for i in a.arange(6):
  s_error_simp_array[i] = abs(Exact_int_s - st.simpson(a.s,0,2,n[i])

在这里，我尝试找出近似值中的错误。问题首先是 Exact_int_s 显然是 -4.5*Si(zoo) + 8.16827746848576，我不知道那是什么意思，而且 simpson 函数总是返回 nan。

我知道有很多文字和代码，但有人知道哪里出了问题吗？

Answer 1

避免回答-4.5*Si(zoo)+ 8.---只需从一个小的正数开始积分，例如:

x = sp.Symbol('x')
print sp.integrate( x * sin(3./x), (x, 0.000001, 2) )

你会得到类似 1.0996940... 的答案

你可以证明这一点，因为 |s(x)| <= x对于小x , 所以区间 [0, epsilon]贡献不了那么多。

顺便说一句 - 您的 simpson 实现似乎已经过检查。