数值分析：梯形公式 、辛普森公式、复化梯形公式、复化辛普森公式
Python实现梯形公式 、辛普森公式、复化梯形公式、复化辛普森公式
数值求积公式概念

梯形公式与辛普森公式

梯形公式与辛普森公式的余项

复化求积公式

复化梯形公式与其余项

复化辛普森公式与其余项

Python实现四种公式.
题目.
Python编写梯形公式、辛普森公式、复化梯形公式、复化辛普森公式

并利用其分别求解sqrt(x) * log(x) 与 sin(x)/x 在(0,1)上的积分。
具体代码实现：
import math
import numpy as np
import matplotlib.pyplot as plt
#待求解数值积分sqrt(x) * log(x)
def f1(x):
    if (float(np.fabs(x))<1e-15) :
        return 0
    y=np.sqrt(x) * np.log(x)
    return y
#待求解数值积分sin(x)/x
def f2(x):
    if (float(np.fabs(x)) < 1e-15):
        return 1
    y=np.sin(x)/x
    return y
#梯形公式 f为待求解积分 a为积分下限 b为积分上限
def TX(f,a,b):
    TX = 0.5 * (b - a) * (f(a) + f(b))
    print("梯形公式计算结果为：TX = ", TX)
#辛普森公式 f为待求解积分 a为积分下限 b为积分上限
def XPS(f,a,b):
    XPS = (b-a)*(f(a)+4*f((a+b)/2)+f(b))/6.0
    print("辛普森公式计算结果为：XPS = ", XPS)
#复化梯形公式 f为待求解积分 a为积分下限 b为积分上限 n为区间等分数
def FHTx(f,a,b,n):
    ti=0.0
    h=(b-a)/n
    ti=f(a)+f(b)
    for k in range(1,int(n)):
        xk=a+k*h
        ti = ti + 2 * f(xk)
    FHTx = ti*h/2
    print("复化梯形公式计算结果为：FHTx = ", FHTx)
#复化辛普森公式 f为待求解积分 a为积分下限 b为积分上限 n为区间等分数
def FHXPs(f,a,b,n):
    si=0.0
    h = (b - a) / (2 * n)
    si=f(a)+f(b)
    for k in range(1,int(n)):
        xk = a + k * 2 * h
        si = si + 2 * f(xk)
    for k in range(int(n)):
        xk = a + (k * 2 + 1) * h
        si = si + 4 * f(xk)
    FHXPs = si*h/3
    print("复化辛普森公式计算结果为：FHXPs = ", FHXPs)

def main():
    a = input("a = ")  # 积分下限
    b = input("b = ")  # 积分上限
    a = float(a)  # 强制转换为float类型
    b = float(b)
    n = input("n = ") #将区间分成为n等份
    n = float(n)
    #TX(f2,a,b) #调用梯形公式求解
    #XPS(f2,a,b) #调用辛普森公式求解
    #FHTx(f2,a,b,n) #调用复化梯形公式求解
    FHXPs(f2,a,b,n) #调用复化辛普森公式求解


if __name__ == '__main__':
    main()

【计算方法数值分析】复化梯形公式、复化辛普森公式和龙贝格数值积分


1、	复化梯形公式


%复化梯形公式
function t=agui_trapz(fname,d2fname,a,b,e)
%fname为被积函数，d2fname为函数fname的二阶导数，a,b分别为下界和上届，e为精度
y=abs(feval(d2fname,a:1e-5:b));
m=max(y);
h=abs(sqrt(12*e(b-a)./m));
n=ceil((b-a)/h)
h=(b-a)/n;
fa=feval(fname,a);
fb=feval(fname,b);
f=feval(fname,a+h:h:b-h+0.001*h);
t=h*(0.5*(fa+fb)+sum(f));




2、复化辛普森公式


%复化辛普森公式
function s=agui_simpson(fname,d4fname,a,b,e)
%fname为被积函数，d4fname为函数fname的四阶导数，a,b分别为下界和上届，e为精度
y=abs(feval(d4fname,a:1e-5:b));
m=max(y);
h=abs((2800*e(b-a)./m).^(0.25));
n=ceil((b-a)/h)
h=(b-a)/n;
fa=feval(fname,a);
fb=feval(fname,b);
s=fa-fb;
x=a;
for i=1:n;
    x=x+h/2;s=s+4*feval(fname,x);
    x=x+h/2;s=s+2*feval(fname,x);
end
s=s*h/6;




3、龙贝格数值积分法


%龙贝格数值积分法
function r=agui_rbg(fname,a,b)
%fname为被积函数，a,b分别为下界和上界
e=1e-6;
i=1;j=1;h=b-a;
T(i,1)=h/2*(feval(fname,a)+feval(fname,b));
T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h))*h/2;
T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);
while abs(T(i+1,i+1)-T(i,i))>e
    i=i+1;
    h=h/2;
    T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h))*h/2;
    for j=1:i
        T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);
    end 
end
T
r=T(i+1,j+1);
end



结果分析



1.由表格比较可得出，无论是复化梯形公式、复化辛普森公式还是龙贝格积分公式，都有着较高的精度，相较之下，龙贝格积分公式精度最高。

2.数值求积的误差除了与被积函数有关之外，还与积分区间的长度(b-a)有关，积分区间越小，则求积公式的截断误差也越小。而在对积分区间作同样的分割的条件下，复合辛普生求积公式比复合梯形公式的计算精度高。

3.在计算速度方面，从表中可看出，复化梯形公式的等分数要比其它两个大得多，且从计算结果上可看出复化辛普生公式也比复化梯形公式的收敛速度快得多。

4.从计算公式复杂性而言，龙贝格数值积分公式的计算量最少。

自适应复化辛普森公式求积算法（C语言实现）
利用复化辛普森公式求积分自适应步骤
基于C语言实现的代码

利用复化辛普森公式求积分自适应步骤
h为步长，a为积分下限，b为积分上限，f为积分函数，n为划分的积分次数，ε为允许误差限（若要求精度为10的负n次方，则误差限为0.5乘以10的负n次方）


基于C语言实现的代码
代码中积分函数f(x)为1/(1+x^2)（在step2中的for循环中可以看到，如需改变函数则改变这里的值）
#include "stdio.h"

void step2(double *h, int *n, double *S, double *S1, double *S2)
{
	int i = 0;
	*S = 0;
	for(i = 0; i < *n; i++)
	{
		*S += ( 2.0/(1 + ((i+0.25)*(*h) - 4) * ((i+0.25)*(*h) - 4) ) - 1/(1 + ((i+0.5)*(*h) - 4) * ((i+0.5)*(*h) - 4) ) + 2.0/(1 + ((i+0.75)*(*h) - 4) * ((i+0.75)*(*h) - 4) ));
	}
	*S2 = 0.5 * (*S1) + (*h)/6.0 * (*S);
}

int step3(double *S1, double *S2)//判断精度是否满足 
{
	if((*S1 - *S2)>0)
	{
		if((*S1 - *S2) < 0.0000005)//0.0000005为10的负6次方精确度 
			return 1;
		else
			return 0;
	}
	else 
	{
		if((*S2 - *S1) < 0.0000005)
			return 1;
		else
			return 0;
	}
}

main()
{
	double h = 8;
	int n = 1;
	int j = 0;
	double S = 0;
	double S1 = 280.0/51.0; 
	double S2 = 0;
	step2(&h, &n, &S, &S1, &S2);
	while(!step3(&S1, &S2))
	{
		h /= 2.0;
		n *= 2;
		S1 = S2;
		step2(&h, &n, &S, &S1, &S2);
	}
	printf("%.10f",S2);
}

复化辛普森公式求二重积分matlab源码
直接拷贝到matlab编辑器，傻瓜式操作。具体算法自行探究，网上都有，小编只提供代码。用的好的请加个关注，篱落～～成殇～～再次先行谢过。
%%%%%%%%%%     2020.6.5        %%%%%%%%%
%%%%%%%%%%复化Simpson公式求二重积分%%%%%%%%%
%%%%%%%%%%     Liu Deping      %%%%%%%%%
clear all;
%%%被积函数及积分上下限导入；
s=input('请输入函数表达式：f = ','s');
f=inline(s);
a = input('请输入积分变量x左边界a的值:');
b = input('请输入积分变量x右边界b的值:');
c = input('请输入积分变量y左边界c的值:');
d = input('请输入积分变量y右边界d的值:');
h1 = input('请输入积分变量x步长h1的值:');
h2 = input('请输入积分变量y步长h2的值:');
m=round((b-a)/h1);
n=round((d-c)/h2);
%%%系数表T，t[i,j]为复化梯形公式的系数；
T=zeros(m+1,n+1);
R=[1,4,1;4,16,4;1,4,1];
for p=0:m/2-1
    for q=0:n/2-1
        for i=1:3
            for j=1:3
                T(i+2*p,j+2*q)=T(i+2*p,j+2*q)+R(i,j);
            end
        end
    end
end
fprintf('系数表：');
T
%%%计算各二维节点函数值，并存放于矩阵F中；
F=zeros(m+1,n+1);
for i=0:m
    for j=0:n
       F(i+1,j+1)=f(a+i*h1,c+j*h2);    
    end
end 
%%%结果输出；
format long;
fprintf('各节点函数值；')
F
fprintf('复化Simpson公式计算结果；')
Snm=(b-a)*(d-c)/(9*m*n)*sum(sum(T.*F))

下面附上一个例题吧，感受一下代码的强大。



只看第(2)问，利用辛普森公式计算积分。(关于第一问的代码可以去关注上一篇文章〈复化梯形公式公式求二重积分matlab源码〉)

步骤:

1.复制代码到matlab编辑器，点击运行;

2.按照提示输入:



3.得到结果(包括系数、节点函数值和积分结果):







如果代码对你有帮助，请加一下关注，如果有问题可以私聊！最后，再次感谢喜欢我的小伙伴。