Maximaで学ぶ逐次代入法（データ点のプロット）

数学の具体的な計算にMaximaを使って、数学もMaximaも同時に学んでしまいましょう。

今回はMaximaを使って数値計算における逐次代入法をみてみたいと思います。

本記事のMaximaの学習事項としては離散データの作成・処理の方法、プロットの方法を学びます。

ただし、図の描画の確認はwxMaxima上でgnuplotを用いた場合で行っています。

インターフェイスなどの動作環境の違いによって適宜変更点があるかもしれません。

逐次代入法
データ点のプロット
平均・分散・標準偏差

逐次代入法

方程式 $x+e^x=0$ は $x_0=0$ として逐次的に数列

$\quad x_{n+1}=−e^{x_n}$

に代入することで、数値的な近似解を求めることができます：

kill(all)$
x[0]:0$
for i: 1 step 1 thru 100 do x[i]:bfloat(-exp(x[i-1]))$
x[100];

から、100回逐次代入して得られる近似解が

$x=-0.5671432904097839$

と求まります。実際

bfloat(x[100]+%e^(x[100]));

から十分な精度で数値的な近似解を与えていることがわかります。

データ点のプロット

この数列が収束していく様子を図示してみたいと思います。

そこで、ステップ数 $n$ と各ステップでの数列の値 $x_n$ の組 $(n,x_n)$ を

データファイル（ファイル名は適当です）result.dat に出力します：

kill(all)$
block(local(x,y), filename: openw("result.dat"),
      n_step: 0, 
      x_plot: 0,
      for i: 0 step 1 thru 100
      do(
         printf(filename, "~f ~f ~%", n_step, x_plot),
         n_step: i+1,
         x_plot: -exp(x_plot)
      ),
      close(filename)
)$

次に、上記データファイルをMaximaで読み込める形に変更します：

xy : read_nested_list ("/保存先ディレクトリのパス/result.dat")$

データを２次元面にプロットします：

plot2d ([discrete, xy],[x, 0, 15],[style,points],
[xlabel, "number of steps"],[ylabel, "value of x_n"]);

以上のようにして作成したのが次の図です。

f:id:pianofisica:20200829111745p:plain

10ステップ目くらいから一定の値に収束しているように見えます。

平均・分散・標準偏差

ここで、過去記事

pianofisica.hatenablog.com

で紹介した、最小二乗法による１次関数でのデータ点の近似を使ってみましょう。

最小二乗法によるフィッティングで得られる１次関数は

$\displaystyle{\ y=Ax+B; \qquad A=\frac{{\rm Cov(n,x)}}{\sigma_n^2}, \qquad B={\rm Avr}_x-\frac{{\rm Cov(n,x)}{\rm Avr}_n}{\sigma_n^2}}$

によって与えられます。ただし公式中の量は以下で与えられます。

２つのデータからなる $N$ 組のデータ $(n,x_n)$ （ $n=1,2,\dots,N$ ）があるとき、

それぞれの成分の平均値を

$\displaystyle{\ {\rm Avr}_n=\frac{\sum_{n=1}^N n}{N}, \qquad {\rm Avr}_x=\frac{\sum_{n=1}^N x_n}{N}}$

で定義します。今の場合、収束していそうな10ステップ目以降のデータを用いることにして

N:100$
Avr_n:sum(xy[i][1],i,11,N)/(N-10);
Avr_x:sum(xy[i][2],i,11,N)/(N-10);

より

$\displaystyle{\ {\rm Avr}_x=54.5, \qquad {\rm Avr}_y=-0.5671272541666755}$

と求まります。

分散は

$\displaystyle{\ {\rm Var}_n=\frac{\sum_{n=1}^N (n-{\rm Avr}_n)^2}{N}=\frac{\sum_{n=1}^N n^2}{N}-{\rm Avr}_n^2 \\ \ {\rm Var}_x=\frac{\sum_{n=1}^N (x_n-{\rm Avr}_x)^2}{N}=\frac{\sum_{n=1}^N x_n^2}{N}-{\rm Avr}_x^2}$

で定義されます。データ例に対しては

Var_n:sum((xy[i][1]-Avr_n)^2,i,11,N)/(N-10);
Var_x:sum((xy[i][2]-Avr_x)^2,i,11,N)/(N-10);

より

$\displaystyle{\ {\rm Var}_x=674.9166666666666, \qquad {\rm Var}_y=8.37474787234071\times10^{-8}}$

と計算されます。

標準偏差は分散の二乗根で定義されます：

$\displaystyle{\ \sigma_n=\sqrt{{\rm Var}_n}, \qquad \sigma_x=\sqrt{{\rm Var}_x}}$

データ例に対しては

sigma_n:sqrt(Var_n);
sigma_x:sqrt(Var_x);

より

$\displaystyle{\ {\rm Var}_x=25.97915831328388, \qquad {\rm Var}_y=2.893915664344887\times10^{-4}}$

と求まります。共分散は

$\displaystyle{\ {\rm Cov}(n,x)=\frac{\sum_{n=1}^N(n-{\rm Avr}_n)(x_n-{\rm Avr}_x)}{N}=\frac{\sum_{n=1}^N n x_n}{N}-{\rm Avr}_n{\rm Avr}_x}$

で定義されます。データ例に対しては

Cov:sum((xy[i][1]-Avr_n)*(xy[i][2]-Avr_x),i,11,N)/(N-10);

より

$\displaystyle{\ {\rm Cov}(X,Y)=-7.194386961787936\times10^{-4}}$

と計算されます。以上より、データ例に対して最小二乗法により

A:Cov/sigma_n^2;
B:Avr_x-Cov*Avr_n/sigma_n^2;

から、次の１次関数を得ます：

$\ \displaystyle{y=-1.065966706278\times10^{-6}x-0.5670691589811834}$

データ点と１次関数を

plot2d([ [discrete, xy], A*x+B], [x, 0,15],
[style, points, lines],[color, blue, red],[legend, "data", "fitting"])$

によってプロットしたのが以下の図です。

f:id:pianofisica:20200829112026p:plain

いまの場合、必ずしも１次近似が必要なわけではありませんが、

数列の収束性を、１次関数の傾きの小ささから実感できるかと思います。

キーワード：Maxima、数値計算、逐次近似、データ作成、プロット