首页 > 解决方案 > Python梯形法作图

问题描述

我正在尝试将梯形方法的结果绘制成图表,我发现要获得可以integrate.trapz(f(x))从 scypy 库中使用的积分的确切值。这是我学习计算数学的第二个月,还不确定图书馆的使用情况。

这是我的代码

from scipy import integrate

def f(x):
    return 1 - x - 4*x**3 + 2*x**5

def trapezoid(a, b, n):
    h = (b - a) / n
    s = (f(a) + f(b))

    i = 1
    while i < n:
        s += 2 * f(a + i * h)
        i += 1
        area = ((h / 2) * s)
    return area

def graphTrapezoid():
    val = []
    err = []
    exact_sum = integrate.trapz(f(b))-integrate.trapz(f(a))
    for i in range (2,100):
        val.append(trapezoid(a,b,i))
        errVal = abs((trapezoid(a,b,i) - exact_sum)/exact_sum)
        err.append(errVal)
    plt.plot(val, err)

a = 1
b = 5
n = 1
graphTrapezoid()   

以下是我得到的错误声明

IndexError                                Traceback (most recent call last)
<ipython-input-134-504c0acc46bf> in <module>
     39 
     40 
---> 41 graphTrapezoid()

<ipython-input-134-504c0acc46bf> in graphTrapezoid()
     29     val = []
     30     err = []
---> 31     exact_sum = integrate.trapz(f(b))-integrate.trapz(f(a))
     32     for i in range (2,100):
     33         val.append(trapezoid(a,b,i))

~\Anaconda3\lib\site-packages\numpy\lib\function_base.py in trapz(y, x, dx, axis)
   4059     slice1 = [slice(None)]*nd
   4060     slice2 = [slice(None)]*nd
-> 4061     slice1[axis] = slice(1, None)
   4062     slice2[axis] = slice(None, -1)
   4063     try:

IndexError: list assignment index out of range

标签: pythonpython-3.x

解决方案


从 scipy.integrate.trapz文档

函数定义为:

scipy.integrate.trapz(y, x=None, dx=1.0, axis=-1)

其中 y 应该是“array_like”。但是,您将两个整数值传递给 trapz 函数。

如果你想在 和 之间进行积分,f也许你可以先定义一个 x 值列表:ba

a = 1
b = 5
step = 0.1 # A smaller step will give you a more precise result
x = []
for i in range(int((b-a)/step + 1)):
    x.append(a + i*step)

然后,对您的函数进行采样:

y = [f(val) for val in x]

最后,您可以整合您的样本:

result = integrate.trapz(y, x, dx=step)

推荐阅读