数值分析第四章节 用Python实现数值积分与数值微分(Matlab补充版)

参考书籍：数值分析 第五版 李庆杨 王能超 易大义编 第4章 数值积分与数值微分
文章声明：如有发现错误，欢迎批评指正

4.1数值积分概论
4.1.1数值积分基本思想
使用某种方法近似替代${\int }_a^{b}f(x)dx=(b-a)f(\xi )$中的$f(\xi )$避免牛顿莱布尼茨公式求积困难。更一般地我们选取区间$[a,b]$的$k$个节点，然后用$f(x_k)$加权平均得到$f(\xi )$近似值如：${\int }_a^{b}f(x)dx\approx \sum _{k=0}^n{A}_{k}f(x_k)$。
4.1.2代数精度的概念
定义1：m次代数精度(m次代数精确度)定义。一般欲使${\int }_a^{b}f(x)\approx \sum _{k=0}^n{A}_{k}f(x_k)$有m次代数精度。要求有$\sum _{k=0}^n{A}_k=b-a,\sum _{k=0}^n{A}_k{x}_k=\frac{1}{2}(b^2-a^2),\dots ,\sum _{k=0}^n{A}_k{x}_k^m=\frac{1}{m+1}(b^{m+1}-a^{m+1})$。

梯形公式

因为只有两个节点$x_0,x_1$，所以只有两个系数$A_0,A_1$。如果满足代数精度为1则有$A_0+A_1=b-a,A_0{x}_0+A_1{x}_1=\frac{1}{2}(b^2-a^2)$。两个方程但四个未知数难以求解所以对$x_0,x_1$作假设。$x_0=a,x_1=b$代入方程组中解得$A_0=A_1=\frac{1}{2}(b-a)$。带入原式有${\int }_a^{b}f(x)dx\approx \frac{b-a}{2}[f(a)+f(b)]$。
$$例子:{\int }_a^b\lambda x{e}^{-\lambda x^{-1}}dx，取a=0，b=1,\lambda =1吧，即{\int }_1^{2}x{e}^{-x^{-1}}dx$$
反正我是积不出来我太菜了。题目来自概率论与数理统计：知$x\sim Exp(\lambda )$求$E(\frac{1}{x})$。我感觉$E(\frac{1}{x})$并不存在可能是我理解错了题目意思。但是这里并不重要，只需使用数值方法求解这个定积分吧。

import math
def function(x):
    return x*pow(math.e,-pow(x,-1))
print("{:.8f}".format((2-1)/2*(function(1)+function(2))))
#0.79047038

function_value=@(x)x*exp(-power(x,-1));
result=(2-1)/2*(function_value(1)+function_value(2));
fprintf('%.8f\n',result);

矩形公式

因为只有一个节点$x_0$，所以只有一个系数$A_0$。如果满足代数精度为1则有$A_0=b-a,A_0{x}_0=\frac{1}{2}(b^2-a^2)$。两个方程且两个未知数直接求解。有$A_0=b-a,x_0=\frac{1}{2}(a+b)$。代入原式有${\int }_a^{b}f(x)dx\approx (b-a)f(\frac{a+b}{2})$。
就同样使用上面那个例子吧。

import math
def function(x):
    return x*pow(math.e,-pow(x,-1))
print("{:.8f}".format((2-1)*function((2-1)/2)))
#0.06766764

function_value=@(x)x*exp(-power(x,-1));
result=(2-1)*function_value((2-1)/2);
fprintf('%.8f\n',result);

4.1.3插值型的求积公式
在区间$[a,b]$选取$n+1$个节点用拉格朗日插值多项式近似替代$f(x)$，即${\int }_a^{b}f(x)dx$变为${\int }_a^b{L}_n(x)$。公式仍然为${\int }_a^{b}f(x)dx\approx \sum _{k=0}^n{A}_{k}f(x_k)$但是$A_k$不再由解方程组给出，为$A_k={\int }_a^b{l}_k(x)dx,k=0,1,\dots ,n$。余项为$R[f]={\int }_a^b[f(x)-L_n(x)]dx={\int }_a^b{R}_n(x)dx$。$R_n(x)=\frac{f^{n+1}(\xi )}{(n+1)!}{\omega }_{n+1}(x),{\omega }_{n+1}(x)=(x-x_0)(x-x_1)\dots (x-x_n)$。
4.1.4求积公式的余项
$R[f]={\int }_a^{b}f(x)dx-\sum _{k=0}^n{A}_{k}f(x_k)=K{f}^{(m+1)}(\eta )$。$K$是不依赖于$f(x)$的待定参数。显然当$f(x)$次数$\le m$时，$f^{m+1}(\eta )=0$所以$R[f]=0$满足m次代数精度。
4.1.5求积公式的收敛性与稳定性
定义2：在求积公式${\int }_a^{b}f(x)dx\approx \sum _{k=0}^n{A}_{k}f(x_k)$中，若$\underset{\begin{array}{c}n\to \infty \\ h\to 0\end{array}}{\lim }\sum _{k=0}^n{A}_k{f}_k(x_k)$。其中 h = max ⁡ 1 ≤ i ≤ n { x i − x i − 1 } h=\max\limits_{1\leq i\leq n}\{x_i-x_{i-1}\} h=1≤i≤nmax​{xi​−xi−1​}，则称求积公式${\int }_a^{b}f(x)dx\approx \sum _{k=0}^n{A}_{k}f(x_k)$是收敛的。
定理2：若求积公式${\int }_a^{b}f(x)dx\approx \sum _{k=0}^n{A}_{k}f(x_k)$中系数$A_k>0(k=0,1,\dots ,n)$，则此求积公式是稳定的。
4.2牛顿-柯特斯公式
4.2.1柯特斯系数与辛普森公式
设将积分区间$[a,b]$划分为n等份，步长$h=\frac{b-a}{n}$，选取等距节点$x_k=a+kh$构造出的插值型求积公式$I_n=(b-a)\sum _{k=0}^n{C}_k^{(n)}f(x_k)$，称为牛顿-柯特斯公式，式中$C_k^{(n)}$称为柯特斯系数。由于这是插值型的求积公式所以有$A_k=(b-a)C_k^{(n)}$并且有$A_k={\int }_a^b{l}_k(x)dx,k=0,1,\dots ,n$(这个之前已经说过)。柯特斯系数的理论求解是可行的因为全部是多项式；但是还是记住考试谁想算啊，记住$\le 4$的应该就好。当$n=1$时为梯形公式有$c_0^{(1)}=c_1^{(1)}=\frac{1}{2}$。

辛普森公式

当$n=2$时为辛普森公式有$c_0^{(2)}=c_2^{(2)}=\frac{1}{6},c_1^{(2)}=\frac{2}{3}$。
就同样使用上面那个例子吧。

import math
def function(x):
    return x*pow(math.e,-pow(x,-1))
print("{:.8f}".format((2-1)*(1/6*function(1)+2/3*function((2-1)/2)+1/6*function(2))))
#0.30860189

function_value=@(x)x*exp(-power(x,-1));
result=(2-1)*(1/6*function_value(1)+2/3*function_value((2-1)/2)+1/6*function_value(2));
fprintf('%.8f\n',result);

记一下$n=3$时的系数

柯特斯公式

当$n=4$时为柯特斯公式有$c_0^{(4)}=c_4^{(4)}=\frac{7}{90},c_1^{(4)}=c_3^{(4}=\frac{16}{45},c_2^4=\frac{2}{15}$。
就同样使用上面那个例子吧。

import math
def function1(x):
    return x*pow(math.e,-pow(x,-1))
def function2():
    cnt=0;lt=[7/90,16/45,2/15,16/45,7/90]
    for i in range(5):
        cnt+=lt[i]*function1(1+i*(2-1)/4)
    return cnt
print("{:.8f}".format(function2()))
#0.77672741

%function1
function output=hanshu1(x)
    output=x*exp(-power(x,-1));
end

%function2
function output=hanshu2()
    output=0;
    lt=[7/90,16/45,2/15,16/45,7/90];
    for i=1:5
        output=output+lt(i)*hanshu1(1+(i-1)*(2-1)/4);
    end
end

%ketesigongshi
fprintf('%.8f\n',hanshu2());

4.2.2偶阶求积公式的代数精度
当阶$n$为欧式时，牛顿-柯特斯公式至少有$n+1$次代数精度。
4.2.3辛普森公式的余项
算就好了，不难算的。我们来推一下一般。前面已经说了K不依赖于$f(x)$，不妨令$f(x)=x^{m+1}$，所以有$K=\frac{{\int }_a^{b}f(x)dx-\sum _{k=0}^n{A}_{k}f(x_k)}{f^{(m+1)}(\eta )}=\frac{{\int }_a^b{x}^{m+1}dx-\sum _{k=0}^n{A}_{k}f(x_k)}{(m+1)!}=\frac{\frac{1}{m+2}(b^{m+2}-a^{m+2})-\sum _{k=0}^n{A}_{k}f(x_k)}{(m+1)!}$。所以有$R[f]=K{f}^{(m+1)}(\eta )$。
4.3复合求积公式
就是一种简单套路。没啥可以拿来说的。
4.3.1复合梯形公式
4.3.2复合辛普森公式
$$例3:f(x)={\int }_0^1\frac{sinx}{x}dx$$

复合梯形公式

from math import sin
def function2(x):
    if x==0:
        return 1
    return sin(x)/x
def function1(n):
    lt=[0+1/n*i for i in range(n+1)]
    cnt=0
    for i in range(n):
        cnt+=1/n*((function2(lt[i])+function2(lt[i+1]))/2)
    return cnt
for i in range(10):
    print("复合梯形公式，{:>2d}个区间，{:>2d}个节点，数值积分为{:>10.8f}".format(i+1,i+2,function1(i+1)))

%func
function result=func(x)
    if x==0
        result=1;
    else
        result=sin(x)/x;
    end
end

%function1
function cnt=function1(n)
    lt=zeros(1,n+1);
    for i=0:n
        lt(i+1)=0+1/n*i;
    end
    cnt=0;
    for i=1:n
        cnt=cnt+1/n*((func(lt(i))+func(lt(i+1)))/2);
    end
end

%fuhetixinggongshi
for i=0:9
    fprintf('复合梯形公式，%2d个区间，%2d个节点，数值积分为%10.8f\n',i+1,i+2,function1(i+1));
end

复合辛普森公式

from math import sin
def function2(x):
    if x==0:
        return 1
    return sin(x)/x
def function1(n):
    lt=[0+1/n*i for i in range(n+1)]
    cnt=0
    for i in range(n):
        cnt+=1/n*(1/6*function2(lt[i])+2/3*function2((lt[i]+lt[i+1])/2)+1/6*function2(lt[i+1]))
    return cnt
for i in range(10):
    print("复合辛普森公式，{:>2d}个区间，{:>2d}个节点，数值积分为{:>10.8f}".format(i+1,i+2,function1(i+1)))

%func
function result=func(x)
    if x==0
        result=1;
    else
        result=sin(x)/x;
    end
end

%function2
function cnt=function2(n)
    lt=zeros(1,n+1);
    for i=0:n
        lt(i+1)=0+1/n*i;
    end
    cnt=0;
    for i=1:n
        cnt=cnt+(1/n)*(1/6*func(lt(i))+2/3*func((lt(i)+lt(i+1))/2)+1/6*func(lt(i+1)));
    end
end

%fuhexinpusengongshi
for i=0:9
    fprintf('复合辛普森公式，%2d个区间，%2d个节点，数值积分为%10.8f\n', i+1,i+2,function2(i+1));
end

就这样吧受不了了。后面好像也不会考。