python代码
 
import numpy as np

a, b = input("积分区间：").split(' ')
n = int(input("子区间个数："))

a = int(a)
b = int(b)
x_i = np.linspace(a, b, n+1) # 区间结点
h = (b-a)/n

# 原函数
def f(x):
    return 4/(1+x**2)

def piece(x_i, i, h): 
    return h/6 * (f(x_i[i-1]) + 4*f((x_i[i-1] + x_i[i])/2) + f(x_i[i]))

S = 0
for i in range(1, n+1):
    S = S + piece(x_i, i, h)

print("复化辛普森积分公式计算结果：" + str(S))
 
输入格式及测试结果
 

利用复化 $Simpson$公式计算积分
 
示例
 0 1 5
 $$I={\int }_0^1\frac{\sin x}{x}dx$$
 
#include <iostream>
using namespace std;

double f(double x){
    if (x == 0) {
        return 1;
    }
    return sin(x)/x;
}

int main(int argc, const char * argv[]) {
    cout<<"依次输入积分区间[a,b],和等分数n"<<endl;
    double a, b;
    int n;
    cin>>a>>b>>n;
    
    double h = (b - a) / n;
    double S = f(a) - f(b);
    double x = a;
    int k;
    for (k = 1; k <= n; ++k) {
        x += h / 2;
        S += 4 * f(x);
        x += h / 2;
        S += 2 * f(x);
    }
    S *= h / 6;
    
    cout<<S<<endl;
    return 0;
}

 
一、实验目的及题目
 
1. 实验目的：
 
(1) 学习用复化辛普森公式及龙贝格方法求解积分并掌握这种方法。
 
(2)了解这些辛普森公式及龙贝格方法的概念，参考课本写出用复化辛普森算法以及龙贝格方法计算目标题目的程序，在matlab 中实现，并用matlab 内置的函数计算出结果，并提出存在的问题。
 
2. 题目：
 
利用复化辛普森公式和龙贝格方法计算下列积分： (1)dx e x ?-5
 
.002
 
(2)dx x x ?20
 
2sin )2sin(cos π
 
二、实验用仪器设备、器材或软件环境 计算机、matlab 软件。
 
三、实验原理、程序框图、程序代码 1.实验原理：
 
根据微积分学基本定理，若被积函数f(x)在区间[a,b]上连续，只要能找到f(x)的一个原函数F(x)，便可利用牛顿-莱布尼茨公式求得积分值。但会经常遇到如下问题：找不到用初等函数，找到了原函数，但因表达式过于复杂而不便计算等等。此时则不能用牛顿-莱布尼茨公式，因此有必要研究如下公式。 1)复化求积公式及原理
 
由于高阶插值的不稳定性，为了提高计算积分的精度，可把积分区间分为若干个
 
小区间，将()I f 写成这些小区间上的积分之和，然后对每一个小区间上的积分应用到辛普森公式，或柯特斯公式，并把每个小区间上的结果累加，所得到的求积公式就称为复化求积公式。
 
辛普森公式的数值积分公式为：
 
?
 
+++-≈
 
b
 
a
 
b f b a f a f a b dx x f )]()2
 
(
 
4)([6
 
)(

1. 算法原理简介
 
步1 将积分区间 [a,b] 分成 n 等分，分点为 xk=a+kh(k=0,1,⋯,n)，其中 h=(b-a)/n。
 步2 记区间 [xk,x(k+1)] 的中点为 x(k+1/2)，在每个小区间上用辛普森公式，
 则得到复化辛普森公式：
 
 
2. 应用实例
 
取 n=10，利用复化辛普森公式计算积分：
 
 
3. 程序代码
 
#include<iostream>
using namespace std;
double fun(double x) 
{	
    return 4 / (1 + x * x);
}

double Simpson(double a, double b, int n) 
{	
    double h = (b - a) / n;
    double s1 = fun(a + h / 2);	
    double s2 = 0;	
    for (int i = 1; i < n; i++) 	
    {  		
        s1 += fun(a + i * h + h / 2);		
        s2 += fun(a + i * h);	
    }	
    return h * (fun(a) + 4 * s1 + 2 * s2 + fun(b)) / 6;
}

int main() 
{	
    int n = 10;  //区间等分数	
    double a = 0, b = 1;  //上下限	
    cout << "被积函数 f = 4 / (1 + x * x)" << endl;	
    cout << "区间等分数为 10" << endl;	
    cout << "积分区间为（0，1）" << endl;	
    cout << "积分值为：" << Simpson(a, b, n) << endl;	
    return 0;
}
 
4. 运行结果