学会了函数，也学会了循环、分支、流程，其实就可以做很多事情了，至少可以拿来解决很多数据相关的问题。比如，求函数的积分。
什么是数值积分？
在数学中，积分是求解函数曲线下面积的重要工具，计算方法可以为解析法和数值法。 简单来说，解析法就是求出原函数，然后来用上下界代入计算，就得到了。 数值法则简单粗暴，就是对付着近似一下，只要近似的精度达到要求就可以。 有很多函数都有原函数，但在实际应用中，有些函数的原函数形式复杂，与其求原函数，不如用数值方法来计算。 再比如可能有的同学像我一样就是不会求原函数，或者懒得算，那我就写代码来算显得厉害等等因素，可以使用数值方法来近似计算积分值。 梯形法就是一种简单而有效的数值积分方法。
梯形法的基本思路
梯形法的核心思想是"化曲为直"，将复杂的曲线用简单的直线来近似：
1.分割区间：将积分区间 [a,b] 分成 n 个等宽的小段2.近似替代：在每个小段上，用梯形面积代替曲边梯形面积3.求和累加：将所有梯形面积相加，得到积分的近似值
数学原理
对于函数 f(x) 在区间 [a,b] 上的定积分，梯形法的公式为：
∫[a,b] f(x)dx ≈ (b-a)/2n * [f(x₀) + 2f(x₁) + 2f(x₂) + ... + 2f(xₙ₋₁) + f(xₙ)]
其中：
•n 是分割的区间数•h = (b-a)/n 是每个小区间的宽度•xᵢ = a + i*h 是第 i 个分割点
梯形面积计算
每个小梯形的面积计算公式为： 第 i 个梯形面积 = [f(xᵢ) + f(xᵢ₊₁)] * h / 2
这其实就是我们熟悉的梯形面积公式：(上底 + 下底) × 高 ÷ 2
Python代码实现解析
让我们通过一个具体的Python代码示例来理解梯形法的实现过程。
被积函数定义
def f(x): """被积函数: f(x) = x * exp(x) - 2""" return x * np.exp(x) - 2
我们选择的被积函数是 f(x) = x * e^x - 2，我们用它来演示梯形法的效果，并与解析解进行对比验证。
该函数的原函数为： F(x) = (x-1) * e^x - 2x + C
因此定积分的解析解为： ∫[a,b] (x * e^x - 2) dx = [(x-1) * e^x - 2x]ₐᵇ = [(b-1) * e^b - 2b] - [(a-1) * e^a - 2a]
梯形法实现步骤
#!/usr/bin/env python3# -*- coding: utf-8 -*-import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationfrom matplotlib.patches import Polygon
# 配置中文字体plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']plt.rcParams['axes.unicode_minus'] = False
def f(x): """被积函数: f(x) = x * exp(x) - 2""" return x * np.exp(x) - 2
def trapezoidal_animate_process(f, a, b, n): """ 参数: f: 被积函数 a: 积分下限 b: 积分上限 n: 分割区间数 """ # 计算分割点和函数值 x_points = np.linspace(a, b, n+1) y_points = f(x_points)
# 计算每一步的梯形面积和累积面积 h = (b - a) / n steps = [] cumulative_area = 0
for i in range(n): # 当前梯形的四个顶点 x_left, x_right = x_points[i], x_points[i+1] y_left, y_right = y_points[i], y_points[i+1]
# 梯形面积计算 trap_area = (y_left + y_right) * h / 2 cumulative_area += trap_area
steps.append({ 'x_left': x_left, 'x_right': x_right, 'y_left': y_left, 'y_right': y_right, 'area': trap_area, 'cumulative': cumulative_area })
# 创建图形 fig, ax = plt.subplots(figsize=(12, 8))
# 绘制完整的函数曲线 x_curve = np.linspace(a-0.5, b+0.5, 1000) y_curve = f(x_curve) ax.plot(x_curve, y_curve, 'b-', linewidth=2, label='$f(x) = x e^x - 2$')
# 填充真实的积分区域 x_fill = np.linspace(a, b, 1000) y_fill = f(x_fill) ax.fill_between(x_fill, 0, y_fill, alpha=0.3, color='lightblue', label='真实积分区域')
# 绘制坐标轴 ax.axhline(0, color='black', linewidth=0.8) ax.axvline(0, color='black', linewidth=0.8) ax.axvline(a, color='green', linestyle='--', linewidth=1.2, label=f'下限 a = {a}') ax.axvline(b, color='green', linestyle='--', linewidth=1.2, label=f'上限 b = {b}')
# 设置图形属性 ax.set_xlim(a-0.5, b+0.5) y_min, y_max = np.min(y_curve), np.max(y_curve) ax.set_ylim(y_min-0.5, y_max+1) ax.set_xlabel('x', fontsize=12) ax.set_ylabel('f(x)', fontsize=12) ax.set_title('梯形法则数值积分动画演示', fontsize=16) ax.legend(fontsize=11) ax.grid(True, alpha=0.3)
# 添加说明文本 info_text = ax.text(0.02, 0.98, '', transform=ax.transAxes, fontsize=11, verticalalignment='top', bbox=dict(boxstyle="round,pad=0.3", facecolor="yellow", alpha=0.7))
# 存储已绘制的梯形 trapezoids = []
def animate(frame): """动画函数""" # 清除之前的梯形（除了第一帧） if frame > 0 and frame <= len(steps): # 移除上一个梯形 if trapezoids: trapezoids[-1].remove() trapezoids.pop()
# 如果是最后一帧，显示所有梯形并停止 if frame >= len(steps): # 显示所有梯形 for step in steps: vertices = [ (step['x_left'], 0), (step['x_left'], step['y_left']), (step['x_right'], step['y_right']), (step['x_right'], 0) ] trap = Polygon(vertices, closed=True, facecolor='red', alpha=0.5, edgecolor='red') ax.add_patch(trap) trapezoids.append(trap)
info_text.set_text(f'梯形法则积分完成\n' f'分割数: {n}\n' f'最终结果: {steps[-1]["cumulative"]:.6f}') return trapezoids + [info_text]
# 获取当前步骤 step = steps[frame]
# 绘制当前梯形 vertices = [ (step['x_left'], 0), (step['x_left'], step['y_left']), (step['x_right'], step['y_right']), (step['x_right'], 0) ] trap = Polygon(vertices, closed=True, facecolor='red', alpha=0.3, edgecolor='red', linewidth=2) ax.add_patch(trap) trapezoids.append(trap)
# 更新信息文本 info_text.set_text(f'第 {frame+1} 个梯形\n' f'区间: [{step["x_left"]:.2f}, {step["x_right"]:.2f}]\n' f'面积: {step["area"]:.6f}\n' f'累积: {step["cumulative"]:.6f}')
return trapezoids + [info_text]
# 创建动画 anim = FuncAnimation(fig, animate, frames=len(steps)+20, interval=1000, repeat=True)
plt.tight_layout() plt.show()
return anim, steps[-1]['cumulative'] if steps else 0
def compare_accuracy(f, a, b, n_values): """ 比较不同分割数下的积分精度
参数: f: 被积函数 a: 积分下限 b: 积分上限 n_values: 分割数列表 """ print(f"\n精度比较 (积分区间 [{a}, {b}]):") print("-" * 50) print(f"{'分割数':<8} {'积分结果':<15} {'误差':<15}") print("-" * 50)
# 计算解析解（用于比较） # ∫(x * e^x - 2)dx = (x-1) * e^x - 2x + C analytical = (b-1) * np.exp(b) - 2*b - ((a-1) * np.exp(a) - 2*a)
for n in n_values: # 数值积分 x = np.linspace(a, b, n+1) y = f(x) h = (b - a) / n numerical = h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])
# 计算误差 error = abs(numerical - analytical)
print(f"{n:<8} {numerical:<15.8f} {error:<15.2e}")
print("-" * 50) print(f"解析解: {analytical:.8f}")
# 主程序if __name__ == "__main__": # 设置参数 a, b = 1, 4 # 积分区间 n = 10 # 分割数 # 创建动画 anim, result = trapezoidal_animate_process(f, a, b, n) print(f"数值积分结果: {result:.8f}")
预览时标签不可点
Close
更多
Name cleared
微信扫一扫赞赏作者
Like the AuthorOther Amount
赞赏后展示我的头像
作品
暂无作品
Like the Author
Other Amount
¥
最低赞赏 ¥0
OK
Back
Other Amount
更多
赞赏金额
¥
最低赞赏 ¥0
1
2
3
4
5
6
7
8
9
0
.
Python语言程序设计 · 目录
Python语言程序设计
上一篇牛顿法求解过程的Python实现及其详解下一篇借助Python来理解自适应辛普森积分：原理、实现与动态可视化
Close
更多
搜索「」网络结果
Close
调整当前正文文字大小
更多
100%