5 サンプル・プログラムの実行例

(2009年7月9日 グラフィックス・ライブラリィを GLSC に乗り換える。)

(2019年1月15日 現象数理学科 Mac ですぐに試せるように書き換えた。)

$A$ は計算しておく。 $\vec F'(U^{(k)})^{-1} \vec F(\vec U^{(k)})$ を 計算するため、 $A_k:=\vec F'(\vec U^{(k)})=A-\diag[U^{(k)}_1,\dots,U^{(k)}_{N-1}]$ とおき、連立1次方程式 $A_k \vec x=\vec F(\vec U^{(k)})$ を解く。 $A_k$ が三重対角行列なので、 Gauss の消去法を使うと非常に効率的に解ける。

入手してコンパイル (GLSC)

% curl -O http://nalab.mind.meiji.ac.jp/~mk/program/misc/Newton-glsc.c
% curl -O http://nalab.mind.meiji.ac.jp/~mk/program/linear/trid-lu.c
% curl -O http://nalab.mind.meiji.ac.jp/~mk/program/linear/trid-lu.h
% cglsc Newton-glsc.c trid-lu.c

実行 (GLSC)

% ./Newton-glsc
H (10位でOK)=10
du=46.0345
du=7.2523
du=0.688334
du=0.00573817
du=4.06013e-07
du=2.50143e-14
min=0, max=11.7958
%	

図 1: Newton 法の近似の様子

(2022/12/10追記) Newton-glsc.c の GLSC3D バージョン Newton-glsc3d.c を作った。

入手してコンパイル(GLSC3D)

curl -O http://nalab.mind.meiji.ac.jp/~mk/program/misc/Newton-glsc3d.c
curl -O http://nalab.mind.meiji.ac.jp/~mk/program/linear/trid-lu.c
curl -O http://nalab.mind.meiji.ac.jp/~mk/program/linear/trid-lu.h
curl -O http://nalab.mind.meiji.ac.jp/~mk/program/misc/ccg2
chmod +x ccg2 ; mv ccg2 ~/bin
ccg2 Newton-glsc3d.c trid-lu.c

実行(GLSC3D)

% ./Newton-glsc
H (10位でOK)=10
du=46.0345
du=7.2523
du=0.688334
du=0.00573817
du=4.06013e-07
du=2.50143e-14
min=0, max=11.7958
%