桂电七院数值分析实验报告四

实验目的

通过实验掌握利用Matlab进行数值积分的操作，深刻理解复化梯形，复化辛普生公式。

实验原理

Matlab中，有内置函数计算积分：
>> z = trapz(x,y)
其中，输入x，y分别为已知数据的自变量和因变量构成的向量，输出为积分值。
>> z = quad(fun,a,b)
这个命令是使用自适应求积的方法计算积分的命令。其中，fun为被积函数，a，b为积分区间。
我们还可以利用复化梯形公式

编写Matlab程序

function s=tixing(a,b,n)  
%输出s为积分的数值解，输入(a,b)为积分区间，n为等分区间的个数.
x=a:(b-a)/n:b;
s=(b-a)/n/2*(f(x(1))+f(x(n+1))); 
for i=2:n
    s=s+(b-a)/n*f(x(i));       
end

注意，首先需要将定积分$$I={\int }_a^{b}f\left(x\right)dx$$
的被积函数$f\left(x\right)$存储为m文件f.m，方可使用以上复化梯形法的程序。

实验内容与步骤

1. 复化辛普生公式的表达式不唯一，但本质相同，首先确定复化辛普生公式的一种表达，编写复化辛普生公式的Matlab的程序。
   复化辛普生如下:
   
2. 利用复化辛普生程序计算$I={\int }_0^1\frac{4}{1+x^2}dx$,与复化梯形公式的情况比较速度。
   在函数运行前加上tic,函数结束之前加上toc，可以查看输出时间,同时设置分为50个小的区间,比较输出时间,输出时间短的速度快。复化梯形公式历时0.000856 秒,复化辛普生公式历时 0.000763 秒,所以复化辛普生公式速度更快。
   3.积分$${\int }_{}^{}\frac{\sin x}{x}dx$$
   的原函数无法用初等函数表达，但可以用积分上限函数表达。结合Matlab复化梯形程序，用描点法绘制其原函数$${\int }_1^x\frac{\sin t}{t}dt$$
   在区间$x\in [2,50]$的图形。
   f.m

function y=f(x)
y=4*(1+x.^2).^(-1);

tixing.m

function s=tixing(a,b,n)  
%输出s为积分的数值解，输入(a,b)为积分区间，n为等分区间的个数.
tic
format long
x=a:(b-a)/n:b;
s=(b-a)/n/2*(f(x(1))+f(x(n+1))); 
for i=2:n
    s=s+(b-a)/n*f(x(i));       
end
toc

Sn.m

function Sn = Sn(a,b,n)      
    format long
    h = (b-a)/n;
    sum1 = 0;
    sum2 = 0;
    for i = 0:n-1
        sum1 = sum1 + f(a+(i+1/2).*h);
    end
    for j = 1:n-1
        sum2 = sum2 + f(a+j.*h);
    end
Sn = h/6*(f(a)+4*sum1+2*sum2+f(b));

tixing1.m

function s=tixing1(a,b,n)  
%输出s为积分的数值解，输入(a,b)为积分区间，n为等分区间的个数.
x=a:(b-a)/n:b;
s=(b-a)/n/2*(g(x(1))+g(x(n+1))); 
for i=2:n
    s=s+(b-a)/n*g(x(i));       
end

question3.m

y=[];
j=0;
for x=2:0.05:50
    j=j+1;
    y(j)=tixing1(1,x,50);
end
x1=2:0.05:50;
plot(x1,y,'-')