for i=1:m, for j=1:m H(i,j)=v(i+j-1); end,end
②考虑某一行(或列)，ai = [hi; hi+1; ￠ ￠ ￠ ; hi+m?1]，就可以用单重循环生成Hankel 矩阵了 function H=myhankel(v)
m=(length(v)+1)/2; % 严格说来还应该判定给定输入向量长度奇偶性 for i=1:m, H(i,:)=v(i:i+m-1); end ③利用现有的hankel() 函数，则 function H=myhankel(v)
m=(length(v)+1)/2; % 严格说来还应该判定给定输入向量长度奇偶性 H=hankel(v(1:m),v(m:end));
第2 章MATLAB 语言程序设计基础9
9 已知Fibonacci 数列由式ak = ak?1+ak?2; k = 3; 4; ￠ ￠ ￠ 可以生成，其中初值为a1 = a2 = 1，
试编写出生成某项Fibonacci 数值的MATLAB 函数，要求 ①函数格式为y=fib(k)，给出k 即能求出第k 项ak 并赋给y 向量； ②编写适当语句，对输入输出变量进行检验，确保函数能正确调用； ③利用递归调用的方式编写此函数。
【求解】假设fib(n) 可以求出Fibonacci 数列的第n 项，所以对n > 3 则可以用k=fib(n ? 1)+fib(n ? 2) 可以求出数列的n + 1 项，这可以使用递归调用的功能，
11
而递归调用的出口为
1。综上，可以编写出M-函数。 function y=fib(n) if round(n)==n & n>=1 if n>=3
y=fib(n-1)+fib(n-2); else, y=1; end else
error('n must be positive integer.') end
例如，n = 10 可以求出相应的项为 >> fib(10) ans = 55
现在需要比较一下递归实现的速度和循环实现的速度 >> tic, fib(20), toc ans = 832040 elapsed_time = 62.0490
>> tic, a=[1 1]; for i=3:30, a(i)=a(i-1)+a(i-2); end, a(30), toc ans =
12
832040 elapsed_time = 0.0100
应该指出，递归的调用方式速度较慢，比循环语句慢很多，所以不是特别需要，解这样问题 没有必要用递归调用的方式。 10 第2 章MATLAB 语言程序设计基础
10 由矩阵理论可知，如果一个矩阵M 可以写成M = A + BCBT , 并且其中A, B, C 为相应
阶数的矩阵，则M 矩阵的逆矩阵可以由下面的算法求出
M?1 =
3
A + BCBT ′?1 = A?1 ? A?1B 3
C?1 + BTA?1B ′?1
BTA?1
试根据上面的算法用MATLAB 语句编写一个函数对矩阵M 进行求逆，并通过一个小例子来
检验该程序，并和直接求逆方法进行精度上的比较。
13
【求解】编写这个函数 function Minv=part_inv(A,B,C)
Minv=inv(A)-inv(A)*B*inv(inv(C)+B'*inv(A)*B)*B'*inv(A); 假设矩阵为
M =
2 664 51 50 36 16 50 77 60 32 36 60 87 48 16 32 48 68 3 775
且已知该矩阵可以分解成
A =
2 664 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 4 3
14
775
; B =
2 664 1 2 3 4 2 3 4 0 3 4 0 0 4 0 0 0 3 775
; C =
2 664 4 0 0 0 0 3 0 0 0 0 2 0 0 0 0 1 3 775
对这个例子。可以
>> M=[51 50 36 16; 50 77 60 32; 36 60 87 48; 16 32 48 68]; iM=inv(M); % 数值逆，直接解法
15