Matlab 的数值积分问题
(1)求和命令sum 调用格式.
如果x 是向量,则sum(x) 给出x 的各个元素的累加和;如果x 是矩阵,则sum(x)是一个元素为x 的每列列和的行向量.
例3.1 调用命令sum 求向量x 的各个元素的累加和。
解:输入
x=[1,2,3,4,5,6,7,8,9,10];
sum(x)
得到
ans=55
例3.2 调用命令sum 求矩阵x 的各列元素的累加和。
解:输入
x=[1,2,3;4,5,6;7,8,9]
x=
1 2 3
4 5 6
7 8 9
sum(x)
得到
ans=12 15 18
2.定积分的概念.
定积分是一个积分和的极限.
例如取x e x f =)(,求定积分⎰10dx e x
的近似值。
积分区间为[0,1],等距划分为20个子区间,
x=linspace(0,1,21);
选取每个子区间的端点,并计算端点处的函数值.
y=exp(x);
取区间的左端点处的函数值乘以区间长度全部加起来.
y1=y(1:20);
s1=sum(y1)/20
s1=1.6757
s1可作为定积分⎰10dx e x 的近似值。
若选取右端点:
y2=y(2:21);
s2=sum(y2)/20
s2=1.7616
s2也可以作为定积分⎰10dx e x 的近似值。
下面我们画出图象.
plot(x,y);hold on
for i=1:20
fill([x(i),x(i+1),x(i+1),x(i),x(i)],[0,0,y(i),y(i),0],'b')
end
如果选取右端点,则可画出图象.
for i=1:20
fill([x(i),x(i+1),x(i+1),x(i),x(i)],[0,0,y(i+1),y(i+1),0],'b')
hold on
end
plot(x,y,'r')
在上边的语句中,for … end 是循环语句,执行语句体内的命令20次,fill 命令可以填充多边形,在本例中,用的是兰色(blue)填充.
可试取50个子区间看一看结果怎样.下面按等分区间计算。
syms k n
s=symsum(exp(k/n)/n,k,1,n);
limit(s,n,inf)
得结果
ans=exp(1)-1
3.计算定积分
例3.6 计算⎰10dx e x .
解:输入命令:
syms x;
int(exp(x),0,1)
得结果
ans=exp(1)-1.
这与我们上面的运算结果是一致的.
⒈ 由给定数据进行梯形求积
假设已经建立起向量T N T N y y y y x x x x ],,,[,],,,[2121 ==,则可用以下语句进行梯形求积:
sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2
MATLAB 提供的trapz()函数也可直接用梯形法求解积分问题,该函数调用格式为 S=trapz(x,y)
[例1-6-17] 试用梯形法求出),0(π∈x 区间内,函数sin(x),cos(x),sin(x/2)的定积分值。
[求解] >> x1=[0:pi/30:pi]'; y=[sin(x1) cos(x1) sin(x1/2)];
x=[x1 x1 x1]; S=sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2
>> S1=trapz(x1,y)
[例1-6-18] 用定步长方法求解积分⎰2
/30)15cos(πdx x 。
[求解] 鉴于求解区域内被积函数有很强的振荡,可先用下述语句绘制被积函数的曲线。
>> x=[0:0.01:3*pi/2,3*pi/2];
y=cos(15*x); plot(x,y)
采用不同的步距,可分别得到积分近似结果。
>> syms x, A=int(cos(15*x),0,3*pi/2) % 求理论值
>> h0=[0.1,0.01,0.001,0.0001,0.00001,0.000001]; v=[]
for h=h0,
x=[0:h:3*pi/2,3*pi/2]; y=cos(15*x); I=trapz(x,y); v=[v;h,I,1/15-I]; end
例12:用不同方法计算12
1x
e dx ⎰ 解 用以下几种方法计算:
1)矩形公式和梯形公式
>> h=0.01; x=1:0.01:2; y=exp(1./x);
z1=sum(y(1:100))*0.01,
z2=sum(y(2:101))*0.01,
z3=trapz(x,y),
z4=trapz(y)*0.01
结果为z1=2.0254,z2=2.0147,z3=2.0201,z4=2.0201. 由此结果可以看出,梯形算法得到的结果z3和 z4为z1和z2的平均值.
2)辛普森公式
法一:先建立M -文件
%M 函数fun1.m
function y=fun1(x)
y=exp(1./x);
>> z5=quad('fun1',1,2)
法二:
>> f=inline('exp(1./x)','x'); z6=quad(f,1,2)
法一中,使用M-文件描述被积函数,法二中使用inline 函数描述被积函数
此外,单变量函数的数值积分还可以采用数值分析中的其他算法进行求解。
基于Simpson 算法,MATLAB 提供了函数quad()、quadl()等函数,具体实现方法可参见MATLAB 帮助系统。
对于双重积分问题的数值解,MATLAB 提供的dblquad()函数可直接求出,由以下示例可见。
[例1-6-19]求解双重积分⎰⎰---+1122
22/)sin(2dxdy y x e x 。
[求解] >> f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');
y=dblquad(f,-2,2,-1,1)
对于更一般的双重积分问题,可通过使用第三方数值积分工具箱实现。
对于三重定积分的数值求解,则可采用MATLAB 提供的triplequad()函数得出,方法类似上例。
[例1-6-20]用数值方法求解三重定积分问题⎰⎰⎰--1000224ππdzdydx xze z y x 。
[求解] >> triplequad(inline('4*x.*z.*exp(-x.*x.*y-z.*z)','x','y','z'), ... 0,2,0,pi,0,pi,1e-7,@quadl)
例13:计算三重积分(sin()cos())x y x z x dxdydz Ω
++⎰⎰⎰.
其中:0x πΩ≤≤,01,11x x ≤≤-≤≤
解:
>> format long; triplequad(inline('x+y*sin(x)+z*cos(x)'),0,pi,0,1,-1,1) ans =
11.869604395451995。