MATLAB求解方程之trapz函数

pulbieup

trapz 是基于梯形法则的离散点积分函数。 调用形式：
! f" O5 Z0 P8 f+ k$ E! W2 z9 PI = trapz(x,y)
( r4 k) s5 {! `' l4 ~/ F* l其中 x 和 y 分别是自变量和对应函数值，以 sin(x) 在 [0,pi] 积分为例：
; W/ j& i% K- W1 f5 H+ k$ m; Wx = linspace(0,pi,1e3);     %生成 [0,pi] 内的一系列离散点
1 T" T/ t! H$ by = sin(x);
I = trapz(x,y)

浮点数误差
由于计算机中都是以二进制形式存储数据，当用十进制数进行计算时，就会存在十进制数二进制数的转换。但是某些十进制数转化为二进制数是一个无限位的小数，这时必须对该无限位小数进行截断计算机才能存储，那么此时就会产生误差，即浮点数误差。
例如十进制的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。
或者在矩阵中寻找数的位置(也相当于是判断两数相等)。
%例
% ]( h4 p" x' g# H>> a = 0.1:0.1:0.5

4 [- ~+ m% q9 k$ xa =

8 K" [1 |+ I' \3 @- H6 G4 p    0.1000    0.2000    0.3000    0.4000    0.5000
' U/ @0 g( D- C5 L1 h
>> find(a==0.3)
. x( U8 G3 ~0 L: W- ~4 p( o
2 r: t7 M$ f  e3 }2 T  F$ \# gans =
  Z+ k( e: {# \! A
4 L8 g2 |+ a0 u2 o   Empty matrix: 1-by-0
- F" F% W* s3 k0 a/ d1 |由于 a 向量中的 0.3 是由 0.1+0.1+0.1 计算得到，在计算过程中就产生了浮点数误差，这也导致在判断 a==0.3 时结果为false，所以 find(a==0.3) 返回一个空矩阵。
! Y6 f0 a( `3 w' o% g在进行数值计算判断两数相等时，最好不要直接判断，而是设立一个容差值，当两个浮点数的差的绝对值小于给定的容差值时，我们就认为这两个浮点数相等。
' _( }2 S9 O( d1 w' A5 n2 y! j比如对于上面的例子：
* r/ O; P/ T) Q( Q>> a=0.1:0.1:0.5
- b& u! ?( s" T; `0 Q3 B
* ^7 s) e4 P& l( @) N+ b7 ]a =
; R+ Z  x1 `2 \8 s' k" o' G
  0.1000   0.2000   0.3000   0.4000   0.5000

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

tol =
" M; ?; Q: ?# N5 O
  5.5511e-15
4 f, A( J9 L2 \6 z8 \* R/ c
>> find(abs(a-0.3)
3 L; V' o0 z6 J" g- U7 Kans =

  3
8 ^+ s1 |$ A* i+ Y
5 w: z3 r: w8 X. S  ]+ |; e
生成一系列有规律名变量
当循环迭代需要把每次迭代结果进行保存时，如果每次迭代的结果是尺寸不同的矩阵，无法用矩阵进行存储，那么可以利用 eval 和 num2str 这两个函数可以生成一系列例如 a1、a2、a3… 变量对结果进行保存(不推荐这种方法，原因是 eval 这个函数有很多缺点)。
' D. I8 }- v* L0 g' N: C5 veval:：将括号内的字符串视为语句并运行。
6 e9 o, j1 H1 M* y8 R* c, G1 F- |num2str： 将数值转换为字符串。
: x  \1 _8 y% f0 ~%例
, M+ M% F) e! g9 X3 jfor ii=1:10
+ p( k0 Q" ~6 N) J/ c    str=['a',num2str(ii),'=1:',num2str(ii)];
% h! K9 p/ p" c    eval(str)
; n! p1 x* H( i+ r8 S0 D6 o# \end
这样可以生成一系列变量 a1、a2…a10 对循环结果进行保存。
* a5 E+ ~' b3 i) L6 S不推荐使用 eval 函数的原因，帮助文档有详细的解释。
4 p& s4 B2 A; ~+ F. D  @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.
3 u3 R; ^- P8 p+ d; W3 C$ eCode within an eval statement can unexpectedly create or assign to a variable already in the current workspace, overwriting existing data.
8 e! S0 c! I6 t9 \Concatenating strings within an eval statement is often difficult to read. Other language constructs can simplify the syntax in your code.
MATLAB 对于这类问题有更好的解决办法，利用元胞数组对结果进行存储。元胞数组是 MATLAB 中的特色数据类型，它的元素可以是任意类型的变量，包括不同尺寸或不同维度的矩阵。 对于上面的例子，利用元胞数组：
for ii=1:10
9 h! m) e3 |7 h  s% j' D    a{ii}=1:ii;
end
8 H' Z& ~* v8 w( o5 A! s) z7 M' A5 i即生成一系列元胞存储循环结果。这样无论是程序的可读性、运行效率还是后续程序对保存结果调用的方便程度，都远胜于 eval 函数。
8 @. e9 r- H0 x$ P) z0 l除此之外，在处理符号变量时如果需要生成一系列有规律名符号变量，例如生成一个多项式：
y = x1+2*x2+3*x3+…+100*x100
eval+num2str 能够实现，但更简便的方法还是利用矩阵：
) ]; w4 l) o' m2 l+ ox=sym('x',[1,100]);
& C5 u: O5 d1 d( [' H3 Fw=(1:100).*x;
. K6 j$ |7 W# _7 _9 ~y=sum(w)
, V" }; ]) N' ]& s4 C
1 Z- o: }& \7 m2 r! c

youOK

MATLAB求解方程之trapz函数