数値解析(プログラム集)

数値解析(Numerical Analysis) 第6学期 2−1−0 の講義で使用するプログラムです.C言語を使用しています. この文章は情報科学科の講義と演習のために 書いたもので,プログラム本体の閲覧は情報科学科からに限定されます.

変数のアドレスとそのポインタを使う簡単なプログラム

variable.c 変数のアドレスとそのポインタを使う簡単なプログラムです. 使い方は, ここ を見てください.

2次方程式の解を求めるプログラム

searchroot1.c 2次方程式の解を求めるプログラム です.
searchroot2.c 丸め誤差による桁落ちを回避することを考えて, searchroot1.c を改良したプログラムです.
使い方は, ここ を見てください.

ニュートン法により方程式の解を求めるプログラム

newton1.c ニュートン法を使って近似解を列挙し, その収束性をみます.
newton2.c, newton3.c 収束判定条件をつけたニュートン法のプログラムです.
使い方は, ここ を見てください. 試してみると, 収束しない場合もあることがわかります. 一般に グラフのゼロ点求める場合 このように gnuplot などを使って さきにグラフを見ておきます. newton3.c では初期値によっては近似解の列が振動して しまい収束しないことがあります. その理由をグラフから考えてください.

連立一次方程式 Ax = b を解く解法

solvematrixeqn.c ガウスの消去法を使って Ax=b の解を求めます. 使い方は, ここを見てください.

三重対角行列のLU分解

tridiagonal.c 係数行列が三重対角行列である場合に Ax =g を解きます.まず A を LU分解し, Ly=g, Ux=y にわけて解きます. 使い方は, ここを見てください.

5/11 の講義において多次元配列とそのアドレスについて説明が間違っていまし た. 次回に訂正しますが, これについて実験(array2.c)してみまし たので興味ある人はコンパイルして実行してみてください.

共役勾配法(the conjugate gradient method)

係数行列 A が正定値な対称行列である場合に Ax=b を共役勾配法を用いて解くのが cgrad.c です. 使い方は以下の通りです.

% cc cgrad.c -lm
% ./a.out
x(0)= 3.000000
x(1)= -1.000000
x(2)= 4.000000
x(3)= 1.000000

と出力されます.

数値積分

trapezoidal.c 台形公式によって数値積分を行うプログラムです.
使い方は, ここ を見てください.

Runge-Kutta 法による常微分方程式の数値解法

diffeqn.c 使い方は以下の通りです.

% cc diffeqn.c -lm
% ./a.out
これで出力結果が output.dat というファイルに書かれます.
% gnuplot
gnuplot> plot "output.dat" using 1:2 with lines
とすれば解曲線のグラフを見ることができます. また以下のようにすればプリンタにグラフを出力できます.
gnuplot>set term postscript
gnuplot>set output "output.ps"
% lpr output.ps
解いている例は, -d^{2}x/dt^{2}=x です.

Foucault(フーコー)の振り子

このプログラム foucault.c では 振り子の振動面が, 地球の自転の影響を受けて北半球においては 時計回りで回転する様子を数値計算し描画します. 常微分方程式を Runge-Kutta 法で解く解法の応用例です. 使い方は以下の通りです.

% cc foucault.c
% ./a.out

これで出力結果が output.dat というファイルに書かれます.

% gnuplot

として

gnuplot> set size 0.721,1.0
gnuplot> plot "output.dat" using 2:3 with lines

とすれば振り子の先端につけた質点の動きを見ることができます. 糸の長さは 9.8 メートル, 重力加速度 9.8 m/s^2, 地球の自転の角速度 omega, 緯度 lambda は, 2omega sin(lambda)=0.05 として計算しています.
プログラムのミスや改良点などがありましたら講義のときか またはメールで私(谷口)に知らせてください. そのような フィードバックを歓迎します.
参考文献(敬称略)
C入門(浦 昭二・原田賢一共編)培風館
現実的なCプログラミング(Steve Oualline著, 岩谷宏訳) ソフトバンク
数値計算の常識(伊理正夫・藤野和建著)共立出版
数値解析入門(山本哲朗著)サイエンス社
数値計算(高橋大輔著)岩波書店
放送大学「応用数学(藤田宏)」第7章「発展系の数値解析」の数値実験プログラム (桂田 祐史著)
力学I(原島鮮著)裳華房

谷口雅治(Masaharu Taniguchi)