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

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()

这里并不能够说完全实现了复化辛普生公式，因为这里面涉及到具体的函数，而我们需要事先知道函数表达式，才能够求出来。
故而这里是以书上的一道例题来写
例题以及书上结果如下：

这里我也并没有按照流程图来写（不过大同小异），而是按照表达式来写
流程图和表达式如下


代码如下：

# coding=gbk;
#因为使用复化辛普生公式会涉及到函数的具体形式
# 所以这里我就暂且令 f(x)=sin(x)/x
import math;
def compute_fx(temp):  #用来计算函数值
    if temp==0:
        return 1.0;
    return math.sin(temp)/temp;

def fuhua_simpson(a1,b1,h1,n1):  #复化辛普生
    S=compute_fx(b1)-compute_fx(a1);
    account=0;
    x=a1; 
    while account<n:
        S+=4*compute_fx(x+h1/2)+2*compute_fx(x);
        x+=h;
        account+=1;
    return S;

list_temp=input("请分别输入积分得上下限以及想要将其几等分：").split(" ");
a=float(list_temp[0]);
b=float(list_temp[1]);
n=float(list_temp[2]);
h=(b-a)/n; #h是步长
result=fuhua_simpson(a, b, h, n);
result=result*h/6;
print(result);

运行结果：

遇事不决，可问春风

#被积函数
def fun(x):
    return x/(4+x*x)
#复合梯形
def tx(a,b,n):
    h=(b-a)/n
    x=a
    s=fun(x)-fun(b)
    for k in range(1,n+1):
        x=x+h
        s=s+2*fun(x)
    result=(h/2)*s
    return result
#复合辛普森
def xps(a,b,n):   
    h=(b-a)/n
    x=a
    s=fun(x)-fun(b)
    for k in range(1,n+1):
        x=x+h/2
        s=s+4*fun(x)
        x=x+h/2
        s=s+2*fun(x)
    result=(h/6)*s
    return result
a=3
b=6
n=9
t=tx(a,b,n)
p=xps(a,b,n)
print(t,p)

0.5620542501164288 0.5619649373692952
代码求的是函数x/(4+x*x)在x=3到6上的积分