Matlab数值积分

实验目的：

1.熟悉Matlab矩阵操作

2.用Matlab实现求积公式：梯形，辛普森，复合梯形，复合辛普森等算法

3.学会数值积分有关应用

实验要求：

1.掌握梯形和辛普森算法

2.掌握复合梯形和复合辛普森算法

3.掌握判断迭代停止的手段

实验内容：

1.矩阵的四则运算、幂运算；矩阵元素的查找、排序、求和、求积差分。

2.用Matlab实现求积公式：梯形，辛普森，复合梯形，复合辛普森等算法

3.学会数值积分有关应用

（补充2，不在实验过程中展示）

实验步骤：

1.矩阵的四则运算，进行加减法的前提是参与运算的两个矩阵或多个矩阵必须具有相同的行数和列数；或者其中有一个或多个矩阵为标量。其加减法定义：，当其中含有标量x时，，且矩阵的加法运算满足“交换律”、“结合律”、“存在零元”、“存在负元”的运算律。矩阵加减法示例：

 　　数与矩阵到乘法、矩阵与矩阵的乘法、矩阵的除法。

 　　C的元素就是用书x乘矩阵A对应的元素而得到，数与矩阵的乘法满足下列运算律：， ， ， 

 　　数乘示例：

 　　A为m*h矩阵，B为h*n矩阵，则两矩阵的乘积为一个矩阵，且且满足“结合律”、“左分配律”、“右分配律”、“单位矩阵的存在性”。矩阵乘法示例：

 　　矩阵的除法是乘法的逆运算，左除“\”和右除“/”。若A和B是标量，A/B和B\A等价。对一般的二维矩阵A和B，当进行A\B运算时，要求A的行数与B的行数相等；当进行A/B运算时，要求A到列数与B的列数相等。矩阵除法示例：

 　　矩阵的幂运算，当矩阵A为方阵是可进行矩阵的幂运算，其定义为：方阵幂运算示例：

 　　矩阵元素的查找，函数find()的作用是进行矩阵元素的查找，它通常与干洗函数和逻辑运算相结合。语法格式：ind=find(X):该函数查找矩阵X中的非零元素，返回这些元素的单下标；[row,col]=find(X,...):该函数查找矩阵X中的非零元素，函数返回这些元素的双下标i和j。

矩阵元素的排序，函数sort()的作用是按升序排序，排序后的矩阵和眼函数矩阵位数一致，语法：B=sort(A):该函数对矩阵A进行升序排序，A可为矩阵或向量；B=sort(A,dim)：当dim=1时矩阵A按列升序，当dim=2时矩阵按行升序；B=sort(...,mode):该函数按mode指定的方式排序（ascend升序，descend降序）。

　　矩阵元素的求和、求积和差分。函数sum()和cumsum()的作用是对矩阵元素求和，函数prod()和cumprod()的作用是对矩阵元素求积，函数diff()的作用计算矩阵的差分示例：

2.求积公式，梯形公式示例：

梯形求积代码：

function t=rctrap(fun,a,b,n)
%复化梯形公式
%n等分
%a,b区间的左右端点
h=(b-a)/n;
t=0;
for i=1:n-1
      t=t+h*feval(fun,i*h+a);
end
t=t+0.5*h*(feval(fun,a+eps)+feval(fun,b))

rctrap

运行部分截图：

辛普森求积代码：

function s=sptrap(fun,a,b,n)
% n，对应的等分点为2n
h=(b-a)/(2*n);
s=0;
for i=1:n
s=s+feval(fun,(2*i-1)*h+a)*4;
end
for i=1:n-1    
s=s+feval(fun,(2*i)*h+a)*2;
end
s=s+feval(fun,a+eps)+feval(fun,b);
s=s*h/3

sptrap

运行部分截图：

根据梯形公式和估计误差公式编写程序计算积分，分别取，用复合梯形公式计算定积分并与精确值比较.然后观察 对计算结果的有效数字和绝对误差的影响.

h=pi/8000;a=0;b=pi/2;x=a:h:b;n=length(x), y=exp(sin(x));
z1=(y(1)+y(n))*h/2; z2=sum(y(2:n-1))*h; z8000=z1+z2,
syms t
f=exp(sin(t)); intf=int(f,t,a,b), Fs=double(intf),
Juewucha8000=abs(z8000-Fs)

Untitled12

运行部分示例：

 　　用复合梯形公式计算定积分，并与精确值比较.然后观察 对计算结果的有效数字和绝对误差的影响.

 　　代码：

function T=rctrap1(fun1,a,b,m)
n=1;h=b-a; T=zeros(1,m+1); x=a; T(1)=h*(feval(fun1,a)+feval(fun1,b))/2;
for i=1:m 
           h=h/2; n=2*n; s=0;
          for k=1:n/2
           x=a+h*(2*k-1); s=s+feval(fun1,x);
end
T(i+1)=T(i)/2+h*s;
end
T=T(1:m);

rctrap1

运行：

用复合辛普森公式计算定积分取n=20001个等距节点，并与精确值比较.然后再取n=13，观察 n对误差的影响.解 由n=2m+1=20001,得m=10000.根据辛普森公式编写并输入下面的程序。代码：

a=0;b=1;m=10000
m=6; h=(b-a)/(2*m); x=a:h:b; 
y=exp((-x.^2)./2)./(sqrt(2*pi));
z1=y(1)+y(2*m+1); z2=2*sum(y(2:2:2*m)); 
z3=4*sum(y(3:2:2*m));
z=(z1+z2+z3)*h/3, syms t,f=exp((-t^2)/2)/(sqrt(2*pi));
intf=int(f,t,a,b), Fs=double(intf), Juewucha=abs(z-Fs)

Untitled13

运行：

 　　复合辛普森数值积分的MATLAB程序comsimpson（fun,a,b,n)。用comsimpson.m和quad.m分别计算定积分取精度为10e-4，并与精确值比较.。代码：

[Q1,FCNT14] = quad(@fun2,0,1,1.e-4,3),
Q2 =comsimpson (@fun2,0,1,10000)
syms x
fi=int(exp( (-x^2)/2)/(sqrt(2*pi)),x,0, 1);
fs=double(fi)
wQ1= double (abs(fi-Q1)), wQ2= double (abs(fi-Q2))

Untitled14

 　　运行：

3.数值分析的应用。利用复合梯形公式求。代码：

function I=tquad(x,y)
n=length(x);m=length(y);
if n~=m
    error(‘向量x,y的长度必须一致‘);
end
h=(x(n)-x(1))/(n-1);
a=[1 2*ones(1,n-2) 1];
I=h/2*sum(a.*y);

tquad

 　　运行：

小结：

在参考PPT和书本的MATLAB编程之后，感觉解决问题使用计算机的方法，就好像使用一种特殊的语言与一堆无机物交流，并通过一些物理的方法使它回应我们。解数学题的编程的关键，就是熟悉数学规律，并使用该编程语言简洁地表达出来。