MATLAB求解方程之trapz函数

pulbieup

7 N' l' o6 q' R7 Btrapz 是基于梯形法则的离散点积分函数。 调用形式：
9 s, D8 M* z( EI = trapz(x,y)
其中 x 和 y 分别是自变量和对应函数值，以 sin(x) 在 [0,pi] 积分为例：
x = linspace(0,pi,1e3);     %生成 [0,pi] 内的一系列离散点
y = sin(x);
I = trapz(x,y)

1 @1 X, h! e) i7 p浮点数误差
1 F0 R1 G! t. j4 h由于计算机中都是以二进制形式存储数据，当用十进制数进行计算时，就会存在十进制数二进制数的转换。但是某些十进制数转化为二进制数是一个无限位的小数，这时必须对该无限位小数进行截断计算机才能存储，那么此时就会产生误差，即浮点数误差。
  r# d3 r% b% p! |8 |1 f, K例如十进制的0.9，在转化为二进制时是无限循环小数0.1110011001100110011...。此时须对该无限位小数进行截断才能保存在内存当中，截断后再转换回十进制时，0.9就变成了0.90000000000000002，这就是浮点数误差的产生过程。
由于浮点数误差的存在，当进行数值计算时就会出现一些不可避免的问题，最常见的就是判断两数相等时得到与预期相反的结果。
: _- _  a8 U7 _% O1 ~$ @  w例：令 a = 0.1+0.2, b = 0.3, 判断 a==b 时，MATLAB 会返回0, 当执行 a-b 时，会发现结果不是精确等于0，而是一个非常小的数5.5511e-17。
, z+ [- A& \5 }4 M" m$ M/ K& E或者在矩阵中寻找数的位置(也相当于是判断两数相等)。
9 x) ~, `* W8 ^/ R$ w6 C- Y%例
>> a = 0.1:0.1:0.5
8 |* ?* N5 @5 r( _. P( W
$ `8 N: {! Y, ~( Va =

    0.1000    0.2000    0.3000    0.4000    0.5000
0 ~& B7 D% k0 B. I$ ?
>> find(a==0.3)
/ m3 z/ w) I. i) L2 n
ans =

   Empty matrix: 1-by-0
由于 a 向量中的 0.3 是由 0.1+0.1+0.1 计算得到，在计算过程中就产生了浮点数误差，这也导致在判断 a==0.3 时结果为false，所以 find(a==0.3) 返回一个空矩阵。
5 J/ D' g: B2 _在进行数值计算判断两数相等时，最好不要直接判断，而是设立一个容差值，当两个浮点数的差的绝对值小于给定的容差值时，我们就认为这两个浮点数相等。
比如对于上面的例子：
>> a=0.1:0.1:0.5
% W; X$ S7 z" Z5 @
' ]1 y% H& O1 |a =
- D5 f2 w- f3 C* B6 i1 h  @, r
  0.1000   0.2000   0.3000   0.4000   0.5000

>> tol=eps(0.3)*10   %设立容差值，一般比这个点的浮点数误差高一到两个数量级即可。eps函数能够求得该点的浮点数误差值。

$ `& R! {, a! ]& Atol =
9 Y) d1 {* O9 K, T' A* a: C
& a& T" y7 `5 H' ~  5.5511e-15
1 b' k( }- R, B3 ]- U) o
>> find(abs(a-0.3)
$ [) s+ ^/ z: S/ Q! |ans =
0 w$ F! o! {1 V
, c. {" S2 r7 n  ]# O# J" C  3
/ K' R1 h; x: k/ `" v% Q; N, c

生成一系列有规律名变量
# q3 w& ?$ S' Z, `4 n当循环迭代需要把每次迭代结果进行保存时，如果每次迭代的结果是尺寸不同的矩阵，无法用矩阵进行存储，那么可以利用 eval 和 num2str 这两个函数可以生成一系列例如 a1、a2、a3… 变量对结果进行保存(不推荐这种方法，原因是 eval 这个函数有很多缺点)。
eval:：将括号内的字符串视为语句并运行。
3 r; o1 u# A! i$ cnum2str： 将数值转换为字符串。
* y. g: A0 {& ^+ [%例
5 B5 W+ ~0 |: V2 a- o- Yfor ii=1:10
    str=['a',num2str(ii),'=1:',num2str(ii)];
$ F4 r0 {& k/ \; _* J$ B% w    eval(str)
end
8 P/ L7 D" [4 z1 f( W这样可以生成一系列变量 a1、a2…a10 对循环结果进行保存。
不推荐使用 eval 函数的原因，帮助文档有详细的解释。
% R+ L5 g8 _  b* i( w1 ZMATLAB® 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.
Code within an eval statement can unexpectedly create or assign to a variable already in the current workspace, overwriting existing data.
( V& y& |" b' q' O& I/ CConcatenating strings within an eval statement is often difficult to read. Other language constructs can simplify the syntax in your code.
2 @" N9 ]: M% {% n, f- W  @9 X& i$ ?6 v& ]MATLAB 对于这类问题有更好的解决办法，利用元胞数组对结果进行存储。元胞数组是 MATLAB 中的特色数据类型，它的元素可以是任意类型的变量，包括不同尺寸或不同维度的矩阵。 对于上面的例子，利用元胞数组：
9 P3 U( ^- }- E6 }for ii=1:10
2 v6 h7 I# E% V( c, V& B    a{ii}=1:ii;
end
3 D3 ?% S9 }: b" F2 M6 n即生成一系列元胞存储循环结果。这样无论是程序的可读性、运行效率还是后续程序对保存结果调用的方便程度，都远胜于 eval 函数。
1 A3 a4 K' ]0 G* J除此之外，在处理符号变量时如果需要生成一系列有规律名符号变量，例如生成一个多项式：
  `6 L( u+ t4 z1 q! G9 sy = x1+2*x2+3*x3+…+100*x100
4 E- L5 ~* V1 {% N) {% Peval+num2str 能够实现，但更简便的方法还是利用矩阵：
7 {! S+ n8 p( M' V* B: I. Mx=sym('x',[1,100]);
w=(1:100).*x;
! N. X. q/ E0 y6 `- z) py=sum(w)

3 d3 E: q$ \7 {7 b' a' A

youOK

MATLAB求解方程之trapz函数