MATLAB求解方程之trapz函数

pulbieup

trapz 是基于梯形法则的离散点积分函数。 调用形式：
I = trapz(x,y)
其中 x 和 y 分别是自变量和对应函数值，以 sin(x) 在 [0,pi] 积分为例：
x = linspace(0,pi,1e3);     %生成 [0,pi] 内的一系列离散点
y = sin(x);
I = trapz(x,y)
. Q+ r* Q2 h1 g" I! N
浮点数误差
# ~3 m) F6 ^; v0 o0 ~/ e, q由于计算机中都是以二进制形式存储数据，当用十进制数进行计算时，就会存在十进制数二进制数的转换。但是某些十进制数转化为二进制数是一个无限位的小数，这时必须对该无限位小数进行截断计算机才能存储，那么此时就会产生误差，即浮点数误差。
例如十进制的0.9，在转化为二进制时是无限循环小数0.1110011001100110011...。此时须对该无限位小数进行截断才能保存在内存当中，截断后再转换回十进制时，0.9就变成了0.90000000000000002，这就是浮点数误差的产生过程。
由于浮点数误差的存在，当进行数值计算时就会出现一些不可避免的问题，最常见的就是判断两数相等时得到与预期相反的结果。
例：令 a = 0.1+0.2, b = 0.3, 判断 a==b 时，MATLAB 会返回0, 当执行 a-b 时，会发现结果不是精确等于0，而是一个非常小的数5.5511e-17。
* T$ L/ F3 w8 ?! R或者在矩阵中寻找数的位置(也相当于是判断两数相等)。
%例
' t2 ~+ v- F. B" ]9 h8 e>> a = 0.1:0.1:0.5
0 }" F% z3 h$ l) |. G
a =

5 a; y5 G4 A$ x. Y    0.1000    0.2000    0.3000    0.4000    0.5000

% a% V9 R( Y: s% s>> find(a==0.3)

5 r. h3 |2 x+ `6 \  w+ o, d) Zans =
1 b( H+ ?  H: f( g/ N" h
& q% m2 m( j7 U0 l   Empty matrix: 1-by-0
由于 a 向量中的 0.3 是由 0.1+0.1+0.1 计算得到，在计算过程中就产生了浮点数误差，这也导致在判断 a==0.3 时结果为false，所以 find(a==0.3) 返回一个空矩阵。
- K+ B, s3 y$ u, }+ s% J在进行数值计算判断两数相等时，最好不要直接判断，而是设立一个容差值，当两个浮点数的差的绝对值小于给定的容差值时，我们就认为这两个浮点数相等。
比如对于上面的例子：
+ l- u' K0 r2 o- g& v8 h>> a=0.1:0.1:0.5

a =
  H# {9 P, s3 C9 j6 M
  0.1000   0.2000   0.3000   0.4000   0.5000

9 N% L4 t: N& g( }>> tol=eps(0.3)*10   %设立容差值，一般比这个点的浮点数误差高一到两个数量级即可。eps函数能够求得该点的浮点数误差值。
+ ~' U4 _' X1 i+ L/ R) n
' p) m; z5 D! T& [- P/ mtol =
" m- g& q  h( b% l
  5.5511e-15

>> find(abs(a-0.3)
ans =

  3


生成一系列有规律名变量
当循环迭代需要把每次迭代结果进行保存时，如果每次迭代的结果是尺寸不同的矩阵，无法用矩阵进行存储，那么可以利用 eval 和 num2str 这两个函数可以生成一系列例如 a1、a2、a3… 变量对结果进行保存(不推荐这种方法，原因是 eval 这个函数有很多缺点)。
eval:：将括号内的字符串视为语句并运行。
num2str： 将数值转换为字符串。
. \( K+ T0 \1 b; C%例
& I/ K( Q8 X  \1 @1 a. Ifor ii=1:10
! X% S: t& d8 r6 C  I    str=['a',num2str(ii),'=1:',num2str(ii)];
    eval(str)
end
' _; n& l3 I" v  a1 p, r这样可以生成一系列变量 a1、a2…a10 对循环结果进行保存。
不推荐使用 eval 函数的原因，帮助文档有详细的解释。
MATLAB® compiles code the first time you run it to enhance peRFormance for future runs. However, because code in an eval statement can change at run time, it is not compiled.
5 H0 M+ `. k( Q3 Q6 NCode within an eval statement can unexpectedly create or assign to a variable already in the current workspace, overwriting existing data.
Concatenating strings within an eval statement is often difficult to read. Other language constructs can simplify the syntax in your code.
+ C6 l+ {+ D# W) j1 _MATLAB 对于这类问题有更好的解决办法，利用元胞数组对结果进行存储。元胞数组是 MATLAB 中的特色数据类型，它的元素可以是任意类型的变量，包括不同尺寸或不同维度的矩阵。 对于上面的例子，利用元胞数组：
( p- W/ @0 L  F# J0 W. |! g3 [0 ffor ii=1:10
' ~& U1 ]$ ?" b    a{ii}=1:ii;
end
即生成一系列元胞存储循环结果。这样无论是程序的可读性、运行效率还是后续程序对保存结果调用的方便程度，都远胜于 eval 函数。
, G+ w) {; \6 I0 o- k0 f3 x除此之外，在处理符号变量时如果需要生成一系列有规律名符号变量，例如生成一个多项式：
y = x1+2*x2+3*x3+…+100*x100
# Q% [. H8 ~* }4 S+ Q9 l& Neval+num2str 能够实现，但更简便的方法还是利用矩阵：
x=sym('x',[1,100]);
w=(1:100).*x;
y=sum(w)

0 B% \/ L/ ]& s' c) w2 _7 V
: E( z( b' j5 ]; c4 _

youOK

MATLAB求解方程之trapz函数