Processing math: 100%

Python数值积分中的辛普森法则建模


辛普森法则是一种非常强大的数值积分工具。为了在需要较少划分的情况下最大限度地提高精度,辛普森法则使用加权因子来计算积分。梯形法则只考虑两个点,𝑥𝑖和𝑥𝑖+1,来估计梯形的面积,但辛普森法则在每一轮中使用多于两个点(或许多条带)(参见下图)。通过加权各种标准来更新不同点的𝑓(𝑥)值,以减少误差。

辛普森1/3法则

这种方法同时计算两条带的面积,因此每次考虑三个不同的“x”值。为了覆盖整个域,条带数n必须是偶数。

对于前两条带,辛普森1/3法则公式如下:

A=13h(f(x0)+4f(x1)+f(x2))

因此,条带配对的总数为:

I=13h((x0)+4f(x1)+f(x2)))+13h(f(x2)+4f(x3)+f(x4))+...

+13h(f(xn4)+4f(xn3)+f(xn2))

+13h(f(xn2)+4f(xn1)+f(xn))

为了简化编程,简化形式可以写成:

I=13h(f(x0)+f(xn))+4(f(x1)+f(x3)+...+f(xn3)+f(xn1))+2(f(x2)+f(x4)+...+f(xn4)+f(xn2))

最后,将各项写成求和形式将得到:

I=13h(f(a)+f(b))+n1i=1,3,54f(xi)+n2i=2,4,62f(xi)

这里,重要的是要记住“n”是偶数。

假设我们想要执行以下积分:

π20xsin(x)dx

现在可以编写编程过程:

导入模块

# Importing module for array from numpy import *

定义需要进行积分的函数

# Defining the function def f(x): return x*sin(x)

输入函数的左右极限(𝑎, 𝑏)以及要将面积划分的梯形数𝑛 (**𝑛必须是偶数**)。

# Defining function domain (i.e. Left and Right limits) a=0 b=pi/2 # Defining number of trapezoids n=4

计算每个梯形的宽度(ℎ = (𝑏 − 𝑎)/𝑛)

# Evaluating the width of trapezoids h=(b-a)/n

创建𝑥数组。由于梯形数为n,所以𝑥的个数为𝑛 + 1,即“𝑥”数组将有𝑛 + 1个元素。这是因为NumPy数组的索引从0开始,而不是从1开始。

# Creating array of x # (Note: For n trapezoids (no. of x's/nodes = n+1)) x=linspace(a,b,n+1)

首先计算12(𝑓(𝑎)+𝑓(𝑏))并将结果赋值给变量“I”。

# First term of I I=0.5*(f(a)+f(b))

运行一个循环,用于计算𝑥 = 1到𝑥 = 𝑛的公式1(因为有𝑛 + 1个节点)。这里需要注意的是,当主循环中的索引j为偶数时,将计算2𝑓(𝑥𝑖),否则将计算4𝑓(𝑥𝑖)项。

# Summing over remaining n-2 trapezoids for j in range(1,n): if j%2==0: I=I+2*f(x[j]) else: I=I+4*f(x[j])

最后,在循环外部,将“I”乘以h/3并将结果赋值给“I”,然后打印结果。

I=(h/3)*I print(f'I = {round(I,6)}')

完整的Python程序如下:

Open Compiler
from numpy import * def f(x): return x*sin(x) # Left and Right limits a=0 b=pi/2 # Number of trapezoids n=4 # Width of trapezoid h=(b-a)/n # array of x # For n trapezoids (no. of nodes=n+1) x=linspace(a,b,n+1) # First term of I I=(f(a)+f(b)) # Summing over remaining n-2 trapezoids for j in range(1,n): if j%2==0: I=I+2*f(x[j]) else: I=I+4*f(x[j]) I=(h/3)*I print(f'I = {round(I,6)}')

输出

对于n=4,程序返回I = 0.999591,这非常接近精确结果。

辛普森3/8法则

它与辛普森1/3法则的区别仅在于考虑了三条带,这意味着每次都包含四个点。前三条带的公式如下:

A=38h(f(x0)+3f(x1)+3f(x2)+f(x3))

因此,所有条带配对的总和为。

I=38[(f(x0)+3f(x1)+3f(x2)+f(x3))+(f(x3)+3f(x4)+3f(x5)+f(x6))+(f(xn3)+3f(xn2+3f(xn1)+f(xn))]

最后将各项写成求和形式,以简化编程。

I=38[(f(a)+f(b))+n2i=1,4,73[f(xi)+f(xi+1)]+n3i=3,6,92f(xi)]

条带的数量n必须是3的倍数。我们将使用两个循环,一个用于第二项,一个用于第三项,因为从上面的公式可以明显看出,第二项从1开始,第三项从3开始。循环计数器将以三递增。

在这里,我们将采用与之前的1/3情况相同的示例:

π20xsin(x)dx

与1/3情况相比,算法的唯一区别在于求解公式2的方式。如果你仔细观察,你会理解到右边第二项的索引从1开始,而第三项的索引从3开始。因此,循环必须分成两部分,如下所示:

# Summing over remaining n-2 trapezoids j=1 while j<n: I=I+3*(f(x[j])+f(x[j+1])) j=j+3 j=3 while j<n: I=I+2*f(x[j]) j=j+3 I=(3*h/8)*I

完整的Python程序如下:

from numpy import * def f(x): return x*sin(x) # Left and Right limits a=0 b=pi/2 # Number of trapezoids n=6 # Width of trapezoid h=(b-a)/n # array of x # For n trapezoids (no. of nodes=n+1) x=linspace(a,b,n+1) # First term of I I=(f(a)+f(b)) # Summing over remaining n-2 trapezoids j=1 while j<n: I=I+3*(f(x[j])+f(x[j+1])) j=j+3 j=3 while j<n: I=I+2*f(x[j]) j=j+3 I=(3*h/8)*I print(f'I = {round(I,6)}')

输出

对于“n=6”,程序返回:

I = 0.999819

Learn Python in-depth with real-world projects through our Python certification course. Enroll and become a certified expert to boost your career.

结论

在本教程中,我们学习了如何在Python中建模辛普森法则以进行数值积分。我们详细阐述并编写了1/3法则和3/8法则的算法。

更新于:2023年10月4日

650 次浏览

开启你的职业生涯

通过完成课程获得认证

开始学习
广告