2021-06-02 matlab自动选取步长的复化梯形求积算法代码及自学笔记

本文链接：https://blog.csdn.net/weixin_51029293/article/details/117488809

matlab自动选取步长的复化梯形求积算法代码及自学笔记

一、思想：用梯形面积近似积分面积，自动划分积分区间（将区间逐次二分），直到误差小于给定的误差限eps

首先，将给定区[a,b]进行n等分，得到n+1个节点 ${X}_k$ （k=0,1,2...n），利用公式（1）计算梯形面积和 ${T}_n$ 近似等于函数积分,需要计算函数在n+1个节点处的函数积分 $f({X}_k)$ ;

将求积区间再次二等分，得到2n+1个节点，新增n个节点，记二分前的某个区间为[ ${X}_k$ , ${X}_{k+1}$ ],（k=0,1,2...n-1），对该区间进行二等分，新增节点为 ${X}_{k+1/2}$ =( ${X}_k$ + ${X}_{k+1}$ )/2,（k=0,1,2...n-1）那么，利用

公式（2）计算梯形面积和 ${T}_{2n}$ 近似等于函数积分，需要计算 $f({X}_{k+1/2})$ 、 $f({X}_k)$ 以及 $f({X}_{k+1})$ ;

${T}_n=\sum_{k=0}^{n-1} \frac{h}{2}[f(X_k)+f(X_{k+1})]=\frac{h}{2}[f(a)+2\sum _{k=1}^{n-1}f(X_k)+f(b)]$ (1）

${T}_{2n}=\sum_{k=0}^{n-1} \frac{h}{4}[f(X_k)+2f(X_{k+1/2})+f(X_{k+1})]=\frac{1}{2}T_n+\sum _{k=0}^{n-1}\frac{h}{2}f(X_{k+1/2})$ (2)

误差记为二分前后两次积分之差，即 $T_{2n}-T_n<\varepsilon$ (3)

且在数学上，该积分极限可表示为 $I-T_{2n}\doteq \frac{1}{3}(T_{2n}-T_n)$ (4)

二、代码实现：下面以函数 $f(x)=\frac{sin(x)}{x}$ 为例，在matlab中进行代码实现，

（1）建立函数文件：

function y=f(x)
if x==0
   y=1,
else
   y=sin(x)/x,
end

（2）写脚本文件

a=input('输入积分下限a=');
b=input('输入积分上限b=');
eps=input('输入误差限eps=');
x=[a,b];
y=f(x);
h=b-a;
T1=h/2*(f(a)+f(b));           %没有进行二等分时梯形面积
while 1                       %表示一直循环下去，直到判断条件成立，跳出循环
      u=h/2;                  %对区间进行二等分
      H=0;                    %给H设定初值0
      x=a+u;                  %x进行二等分后从a以相应的步长u向右移动直到b 
      while x<b               %该内循环用以计算二分后节点对应的函数值，见（2）
      H=H+f(x);
      x=x+h;                  %此处令x=x+h是为了让x从a以相应的步长u向右移动直到b后跳出循环
      end
       T2=1/2*(T1+h*H),       %（2）
        if abs(T2-T1)<eps     %判断
        I=T2+(T2-T1)/3;       %数学上的近似计算式
        disp(I),              %表示显示I的值
        break;                %跳出循环
        end
     h=u,T1=T2                %此处令h=u，是为了上面的u=h,也就是再次二分区间，此时将T2的值赋给T1，以此迭代
end

三、附赠matlab求积分的函数一枚

int函数 int（f(x),x,a,b）其中f(x)为被积函数，x为积分对象，a,b分别为积分下限和积分上限，使用该函数需进行强制转化

以f(x)=sin(x)/x在[0,1]积分为例

vpa(int(sin(x)/x,x,0,1)        %该结果即为函数sin(x)/x在[0,1]上的积分结果

四、总结：（1）求和问题用while do循环嵌套if break 判断跳出循环比较好用

（2）复杂问题分解简单化，将一个复杂问题拆成多个简单的小问题

（3）今天也要加油呀！

注：博文仅为本人学习所记，如有错误，欢迎各位大佬批评指正。