用Python实现数值积分

1. 引言

本文主要用于对比使用Python来实现数学中积分的几种计算方式,并和真值进行对比,加深大家对积分运算实现方式的理解。

闲话少说,我们直接开始吧。 :)

2. 举个栗子

为了加深大家的印象,首先我们来看个例子:
对下述三次函数
f ( x ) = x 3 − 4 x 2 + 4 x + 2 f(x)=x^3-4x^2+4x+2 f(x)=x34x2+4x+2
进行积分运算:
I = ∫ a b f ( x ) d x = ∫ a b ( x 3 − 4 x 2 + 4 x + 2 ) d x   I = \int_a^bf(x)dx= \int_a^b(x^3-4x^2+4x+2)dx\, I=abf(x)dx=ab(x34x2+4x+2)dx
图示如下:
在这里插入图片描述
如果 a = 1 / 2 a=1/2 a=1/2 , b = 5 / 2 b=5/2 b=5/2,则上述积分的值为:

I = ( x 4 4 − 4 x 3 3 + 2 x 2 + 2 x ) ∣ a b = 61 12 ≈ 5.0833 I=(\frac{x^4}{4}-\frac{4x^3}{3}+2x^2+2x)|_{a}^{b}=\frac{61}{12} \approx 5.0833 I=(4x434x3+2x2+2x)ab=12615.0833

3. 矩形计算面积

我们知道,在数学中,积分运算表示上述曲线和x轴围成的封闭区域的面积,为此,我们在数值预算中,来近似计算上述区域的面积,最直观的想法就是拆成一个个小的矩形来计算对应的面积。

3.1 左侧边长计算面积

为了计算每个小矩形的面积,设计到边长高的选择,这里我么以左侧函数取值作为对应矩形的高来计算相应的小矩形的面积,图示如下:
在这里插入图片描述

对应的代码如下:

import numpy as np
x = np.linspace(0, 3, 1001)
f = lambda x: x**3 - 4*x**2 + 4*x + 2
a = 0.5
b = 2.5
Ax = np.linspace(a, b, 101)
Ay = f(Ax)
def defInt_left(f, a, b, N):
    # left-hand point
    result = 0; FX = []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += f(a)*dx
        FX += [f(a)]
        Xn += [a]
        a += dx
    return result, FX, Xn, dx
N = 4
I_left, FX, Xn, dx = defInt_left(f, a, b, N)
print(I_left)

上述代码中,我们将横坐标拆分为4小份,也就是拆分成4个小矩形,然后使用函数左侧的点坐位小矩形的高,上述代码的运行结果如下:

5.25

3.2 右侧边长计算面积

这里和上述原理类似,只不过每个小矩形的高采用右侧边长函数取值来近似计算,图例如下:
在这里插入图片描述
样例代码如下:

def defInt_right(f, a, b, N):
    # right-hand point
    result = 0; FX = []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += f(a + dx)*dx
        FX += [f(a + dx)]
        Xn += [a]
        a += dx
    return result, FX, Xn, dx

N = 4
I_right, FX, Xn, dx = defInt_right(f, a, b, N)
print(I_right)

运行结果如下:

5.0

3.3 中值边长计算面积

看了上述两种近似计算方式,有同学就说有取左侧点算出来面积大的,有取右侧点算出来面积小的,那干脆折中一下,我们来以中值坐位矩形的高来计算对应的面积。图例如下:
在这里插入图片描述
代码实现如下:

def defInt_middle(f, a, b, N):
    # middle point
    result = 0; FX = []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += f(a + dx/2)*dx
        FX += [f(a + dx/2)]
        Xn += [a]
        a += dx
    return result, FX, Xn, dx

N = 4
I_mid, FX, Xn, dx = defInt_middle(f, a, b, N)
print(I_mid)

运行结果如下:

5.0625

4. 梯形计算面积

读到这里的同学可能会思考,既然可以将封闭区域划分成一个个的小矩形,那当然也可以将其划分成梯形来近似计算相应的面积,图例如下:
在这里插入图片描述
样例代码如下:

def defInt_trapezoid(f, a, b, N):
    # trapezoidal rule
    result = 0; FXa, FXb = [], []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += (f(a) + f(a + dx))*dx/2
        FXa += [f(a)]; FXb += [f(a + dx)]
        Xn += [a]
        a += dx
    return result, FXa, FXb, Xn, dx

N = 4
I_trap, FXa, FXb, Xn, dx = defInt_trapezoid(f, a, b, N)
print(I_trap)

运行结果如下:

5.125

5. 真值比对

最后,我们来针对不同的N来讲封闭区域划分成对应的小份,分别针对性的计算上述四种方式的积分值,样例代码如下:

Nx = range(1, 11)
I1, I2, I3, I4 = [], [], [], []
for Ni in Nx:
    i1, *_ = defInt_left(f, a, b, Ni); I1 += [i1];
    i2, *_ = defInt_right(f, a, b, Ni); I2 += [i2];
    i3, *_ = defInt_middle(f, a, b, Ni); I3 += [i3];
    i4, *_ = defInt_trapezoid(f, a, b, Ni); I4 += [i4];

最后将其与真值进行对比,如下:

在这里插入图片描述
可以看出,随着划分区域的增多,采用梯形计算面积方式最逼近真值。

6. 总结

本文重点介绍了使用不同面积划分方法来近似计算积分取值的原理和相应的代码实现,其中采用梯形计算面积的方式随着划分子区域数目的增加最接近真值,推荐大家使用。

您学废了吗?

在这里插入图片描述
关注公众号《AI算法之道》,获取更多AI算法资讯。

  • 3
    点赞
  • 18
    收藏
  • 打赏
    打赏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:大白 设计师:CSDN官方博客 返回首页
评论 1

打赏作者

赵卓不凡

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值