数值分析：梯形公式 、辛普森公式、复化梯形公式、复化辛普森公式

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.从计算公式复杂性而言，龙贝格数值积分公式的计算量最少。