7.5 INTLAB を使わないで試せること

プログラムの入力・編集はこんな感じ

>> edit prog01.m

Kulisch-Miranker [19] の例である。

prog01.m

format long
A=[64919121-159018721;41869520.5-102558961
b=[1;0]
x=A\b;

>> prog01
>> x

>> A*x

Rupm の例 (Loh-Walster [20])

prog02.m

format long;
a=77617; b=33096;
f=(333.75-a^2)*b^6+a^2*(11*a^2*b^2-121*b^4-2)+5.5*b^8+a/(2*b)

次は Foster [21] の例。

http://www.math.sjsu.edu/~foster/publications.htmlからgfpp.mを入手して ²、 ~/Documents/MATLAB に保存する。

N=100;
A=gfpp(N,5);
y=rand(N,1); b=A*y;
x1=A\b;
norm(x1-y)/norm(y)

この結果は “ans = 2.734469866454044e+11” のようになった。とんでもない誤差である。

(epsは浮動小数点数の相対精度。 “This MATLAB function returns the distance from 1.0 to the next larger double-precision number, that is, eps = 2^-52.” と書いてある。

   eps$$=\min\left\{x\in\mathbb{F}\mid x>1\right\}-1=2^{-52} =2.220446049250313080847263336181640625\times10^{-16}$$

ということだよね。 MATLAB の日本語マニュアルでは 「関数 eps は、1.0 から次の最大倍精度数値の距離を返します。つまり、eps = 2^(-52) です。」 となっている。ひどい誤訳だ。)

x2=gmres(A,b,N,N*eps,N);

この結果は “gmres は、相対残差 2e-14 をもつ解に 反復 45 で収束しました。” のようになった。

norm(x2-y)/norm(y)

この結果は “ans = 6.406927643920526e-13” となった。