12.2 自分で関数を作ろう

MATLAB, Octave にはあらかじめ豊富な関数が揃っているが、自分でも作れる。

関数は、関数名 + ``.m'' という 名前のファイルにその定義を書くことで定義できる。

次に掲げるのは前回配ったプログラム (に少し注釈を加えたもの) である。

cg1.m

   1 % cg1.m --- CG 法で A x = b を解くための関数 cg1()
   2 %  2003/5/15 作成, 5/21 修正
   3 %
   4 % 【使い方】
   5 %    x = cg1(A, b, 初期ベクトル, 相対残差)
   6 
   7 % 【例】
   8 %  まず問題を用意
   9 %   n=2;a=[2 1;1 3];exact=ones(n,1);b=a*exact
  10 %  初期ベクトルを 0 ベクトル、相対残差を 10^{-12} にして実行
  11 %   cg1(a,b,zeros(n,1),1e-12)
  12 %
  13 
  14 function x = cg1(a,b,x0,epsilon)
  15   x = x0;
  16   r = b - a * x;
  17   p = r;
  18   eps_b = epsilon * norm(b);
  19   k = 0;
  20   while norm(r) > eps_b
  21     ap = a * p;
  22     pap = ap' * p;
  23     alpha = p'*r/pap;
  24     x = x + alpha*p;
  25     r = r - alpha*ap;
  26     beta = - (ap'*r)/pap;
  27     p = r + beta * p;
  28     k=k+1;  % 最初 k++ としていましたが、それは MATLAB ではダメです
  29     % 以下の3行は中間結果の表示用 (不要ならば削除すればよい)
  30     k
  31     x
  32     r
  33   end

• 縦ベクトル x, y の内積は y'*x ( $\because x\cdot y=y^T x$ )
• ベクトル x の Euclid ノルム は norm(x) (ちなみに $p$ ノルムは norm(x,p) で、 $\max$ ノルムは norm(x,inf)