x - 牛と戯れた日々を想って

GAMBAS 始めました。
Linux で動かす BASIC プログラム
第 4 回　数値計算
by HOLY
Ver.0（2014.10.03〜2015.06.27）
はじめに
ずいぶん前からコンピュータは数値計算や研究活動の道具と
して使われてきました。現在、コンピュータは市場のニーズに基づ
いて新しい利用法を以って発展していると言えます。
つまり、情報の伝達と通信です。
端末はかつて建物の１フロアを占めるほどの大きさから、ノート
型になり、タブレットや携帯型になっています。現在、生活の場面
では、コンピュータが数値計算の道具であることが忘れられよう
としています。
８０年代頃のパソコンの入門書にはたいてい BASIC のサンプ
ルプログラムが掲載されておりました。ハナタレ小僧の頃、それら
のリストの中ではゲームばかり目がいってたんですが、その中で
もキラリと光る分野がありました。そう、数値計算です。今にしてみ
れば、それら数値計算のサンプルは重要なことを書いてたなあと
想います。そりゃあ当時、電子計算機の本ってのは、その筋のプロ
たちが書いてたんですからねぇ。
この章ではコンピュータの基本的な役割にたち帰って、数値解
法を紹介したいと想います。
(仮称)十進 BASIC と比較
GAMBAS の結果を比較するために(仮称)十進 BASIC のリス
トと結果も記述することにします。(仮称)十進 BASIC は、JIS 規格
に準拠しており、標準規定の BASIC 言語として比較するのにちょ
うどよい位置づけです。
(仮称)十進 BASIC も自由にダウンロードして使用することので
きる BASIC 言語のひとつです。しかもれっきとした日本製で、初
学者や初心者に限らず、ベテランまで普段から電卓代わりに使え
るツールとして国内では知られた由緒正しい BASIC です。
プラットフォームも、Windows、MaxOS、Linux のように主な OS
を網羅してるので、複数端末をもってるとソースを使いまわせるの
でありがたいです。
両方のプログラムを比較してみると、 (仮称)十進 BASIC の特
徴は
（１） 変数の型の宣言がいらない。
（GAMBAS は変数宣言が必要です）
（２） (仮称)十進 BASIC では while〜end 文はない。
（Loop 文は、For〜Next と Do〜until がある）
（３） 行番号なしの構造化 BASIC 様の書き方でも、行番号を使う
方法にも対応。（昔から使ってるユーザへのアフターケア的仕様）
（４） 十数年間、ずっと名前に”(仮称)”がついたままである。
（以下、(仮称)十進 BASIC を JIS-BASIC と表記します）
平方根
加減乗除の演算のみを用いて正の数の平方根の近似値を求
める方法を考えます。
縦の長さが a、横の長さが 1 の長方形と面積の等しい正方形
の 1 辺の長さが √ a です。そこではじめの長方形を面積を変えな
いようにして縦と横の長さが近づくように変形していけばこの長
方形の辺の長さは √ a に近づくはずです。
変形して得られる長方形の縦の長さをもとの長方形の縦、横
の長さの平均にとることにします。すなわち、第 k 回目に得られる
長方形の縦横の長さをそれぞれ x(k), y(k)とするとき、
x( k +1)=
x( k) + y( k )
2
とします。すると、この変形で面積を変えないという条件から、横
の長さ y(k+1)は、
y ( k +1)=
a
x ( k +1)
として定められます。これらの関係式から x(k)のみの式を導くと、
1
a
x(0) =a, x( k +1)= ×( x ( k ) +
)
2
x(k)
となります。次のプログラムは、a を入力すると、この数列{x(k)}の
第 8 項までを求めます。例は a=2 で、2 の平方根 1.41421356・・・
に収束していく様子です。
(JIS-BASIC)
10 INPUT a
20 LET x=a
30 FOR n=1 TO 8
40 LET x=(x+a/x)/2
50 PRINT n,x
60 NEXT n
70 END
(Result)
?2
1
2
3
4
5
6
7
8
1.5
1.41666666666667
1.41421568627451
1.41421356237469
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
(GAMBAS)
' Gambas module file
Sub Main()
Dim a As Float
Dim x As Float
Dim n As Integer
Print " root(a): a= ? ";
Input a
x=a
For n = 1 To 8
x = (x + a / x) / 2
Print n, x
Next
End
(Result)
root(a): a= ? 2
1
1.5
2
1.41666666666667
3
1.41421568627451
4
5
6
7
8
1.41421356237469
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
たとえば SQR(36)の計算結果は正確に 6 になります。また角の
大きさの単位が度である場合には SIN(30)の計算結果は 0.5 に
なります。しかし計算結果の真の値が 16 桁の十進数で表現でき
ないときその結果は近似値になります。たとえば 1/3 の計算結果
には誤差が含まれます。
GAMBAS の数値変数は 15 桁の精度を持ちます。有効数字の
桁数が 15 桁を越える数値を代入すると有効数字が 15 桁になる
ように数値が丸められます。
一方、JIS-BASIC の数値式の計算結果は 16 桁を越える精度
を持ちます。たとえば 1/3 の計算結果は
0.333333333333333333333333333 となって有効数字が 27 桁あ
ります。ただし PRINT 文は有効数字が 15 桁になるように計算結
果を丸めて表示します。
数値式の桁数が数値変数の桁数より大きいのは日本工業規
格(JIS)の規定によるものですが、そこが GAMBAS の計算結果と
は異なる現象を引き起こします。たとえば次のプログラムでは
JIS-BASIC だと 70 行目の実行結果は 0 ではありません。これは
JIS-BASIC での数値式と数値変数の丸める桁数に相違がある
ためです。
（JIS-BASIC）
10 Print " SQRT(36)="; Sqr(36)
20 PRINT " sin(30)="; SIN(PI / 6)
30 ! 数値の丸め
40 LET A=1/3
50 Print A
60 Print 1 / 3
70 PRINT 1/3-A
80 END
(Result)
SQRT(36)= 6
sin(30)= .5
.333333333333333
.333333333333333
3.33333333333333E-16
(GAMBAS)
' Gambas module file
Sub Main()
Dim A As Float
Print " SQRT(36)="; Sqr(36)
Print " sin(30)="; Sin(Pi / 6)
' 数値の丸め
A=1/3
Print A
Print 1 / 3
Print 1 / 3 - A
End
(Result)
SQRT(36)=6
sin(30)=0.5
0.33333333333333
0.33333333333333
0
桁落ち
（JIS-BASIC）
10 LET A=1/3
20 PRINT A-0.33333333333333
30 END
(Result)
.000000000000003
(GAMBAS)
' Gambas module file
Sub Main()
Dim A As Float
A=1/3
Print A - 0.33333333333333
End
(Result)
3.3306690738755E-15
JIS-BASIC を実行すると 20 行では A の値が
0.333333333333333
になっているのですから実行結果は
.000000000000003
になります。けれども、この結果は 3.33333333333333E-15 と書か
れてもいいわけで、この結果は有効数字が 1 桁になってしまった
といえます。
一方、GAMBAS では、計算結果は・・・33 と無限の循環少数に
なるはずですが、e-17 からおかしな数値になっています。
このように計算結果が 15 桁では正確に表現できない数値の
引き算を行うと、正確に表現できなくなる現象を「桁落ち」と呼ん
でいます。
関数 y=f(x)の x=a における微分係数 f'(a)は増分 h を限りなく
0 に近づけたときの平均変化率 f ( x +hh)−f ( x ) の極限で定義されま
す。次のプログラムは f ( x )= √ x であるときに x=2 における平均変
化率が h の値を 10E-1 10E-2 10E-3 10E-4 ... 10E-15 と変化さ
せたときにどう変化するか調べようとしたものです。h を 0 に近づ
けたときの極限は理論的な計算では 2 1 2 (≒0.353553390593274)で
√
すが実行結果をみると h が 0 に近づくとかえって正しい値から離
れていってしまっています。これは上に述べた減算に伴う桁落ち
によるものです。
（JIS-BASIC）
10 DEF f(x)=SQR(x)
20 FOR i=1 TO 15
30 LET h=10^(-i)
40 PRINT h, (f(2+h)-f(2))/h
50 NEXT i
60 END
(Result)
.1
.01
.001
.349241122458488
.35311255026876
.3535092074646
.0001
.353548971286
.00001
.35355294866
.000001
.3535533464
.0000001
.353553386
.00000001
.3535534
.000000001
.3535534
.0000000001
.353554
.00000000001
.35356
.000000000001
.3536
.0000000000001
.354
.00000000000001
.36
.000000000000001
0
(GAMBAS)
' Gambas module file
Sub Main()
Dim i As Integer
Dim h As Float
For i = 1 To 15
h = 10 ^ (- i)
Print h, (Sqr(2 + h) - Sqr(2)) / h
Next
End
(Result)
0.1 0.34924112245849
0.01 0.35311255026873
0.001 0.35350920746446
0.0001 0.35354897128581
1E-5 0.35355294865091
1E-6 0.35355334637721
1E-7 0.35355338434684
1E-8 0.35355338656728
1E-9 0.35355340877175
1E-10 0.35355274263793
1E-11 0.35353941996163
1E-12 0.35349501104065
1E-13 0.3530509218308
1E-14 0.35527136788005
1E-15 0.22204460492503
円周率 π
直径 1 の円に内接する正 n 角形の辺の長さの和は
ですから、正 2k 角形の辺の長さの和を ak とすると、
となります。
半角の公式から
cos
√
180 °
1
180 °
= ⋅(1+cos k −1 )
k
2
2
2
180 °
n
180 °
ak =2k⋅sin
k
2
となるので、
n⋅sin
cos90 °, cos45° , cos
が順に計算できます。
sin
45 °
180 °
,⋯ , cos k , ⋯
2
2
√
180 °
180 °
= 1−cos2 k
k
2
2
の関係を利用すれば
ak の計算ができるのですが、実際にはうまくいきません。k が大き
くなると、 180k ° が 0 に近くなるので、 cos 180k ° は 1 に近くなりま
2
2
2
す。そのため、
1−cos2
180 °
2k
の計算で桁落ちが発生します。
そこで、2 倍角の公式を変形して得られる関係式
sin
θ
sinθ
=
2
θ
2⋅cos
2
を用いて、
cos90 °, cos45° , cos
45 °
180 °
,⋯ , cos k , ⋯
2
2
を計算することにします。
cos 90° =1、sin 90° = 0 を出発点としてこれらの値を順に計算
し 2 k と 2k sin 180k ° を出力するプログラムを作ると次のようになり
2
ます。
（JIS-BASIC）
10 LET c=0
20 LET s=1
30 FOR k=2 TO 30
40 LET c=SQR((1+c)/2)
50 LET s=s/c/2
60 PRINT 2^k,2^k*s,pi-2^k*s
70 NEXT k
80 END
(Result)
4
8
16
32
64
128
256
512
2.82842712474619
3.06146745892071
3.12144515225805
3.13654849054594
3.14033115695475
3.14127725093277
3.14151380114429
3.14157294036708
.313165528843605
8.01251946690812E-2
2.01475013317452E-2
5.04416304385724E-3
1.26149663504124E-3
3.15402657018838E-4
7.8852445498838E-5
1.9713222709718E-5
1024
2048
4096
8192
16384
32768
65536
131072
262144
524288
1048576
2097152
4194304
8388608
16777216
33554432
67108864
134217728
268435456
536870912
1073741824
3.14158772527715
3.14159142151119
3.14159234557011
3.14159257658486
3.14159263433856
3.14159264877698
3.14159265238659
3.14159265328899
3.14159265351459
3.14159265357099
3.14159265358509
3.14159265358862
3.1415926535895
3.14159265358972
3.14159265358977
3.14159265358978
3.14159265358979
3.14159265358979
3.1415926535898
3.1415926535898
3.1415926535898
4.928312639958E-6
1.232078604758E-6
3.08019685846E-7
7.700493359E-8
1.9251235286E-8
4.812812349E-9
1.203207434E-9
3.00802928E-10
7.5201802E-11
1.8798899E-11
4.700794E-12
1.177579E-12
2.96775E-13
7.4477E-14
2.4145E-14
1.0724E-14
7.368E-15
7.368E-15
-6.054E-15
-6.054E-15
-1.142281792E-14
(GAMBAS)
' Gambas module file
Sub Main()
Dim c As Float
Dim s As Float
Dim k As Integer
c=0
s=1
For k = 2 To 30
c = Sqr((1 + c) / 2)
s=s/c/2
Print 2 ^ k, 2 ^ k * s, Pi - 2 ^ k * s
Next
End
(Result)
4
8
16
32
64
128
256
512
1024
2048
4096
8192
16384
32768
2.82842712474619
3.06146745892072
3.12144515225805
3.13654849054594
3.14033115695475
3.14127725093277
3.1415138011443
3.14157294036709
3.14158772527716
3.1415914215112
3.14159234557012
3.14159257658487
3.14159263433856
3.14159264877699
0.3131655288436
0.08012519466908
0.02014750133174
0.00504416304385
0.00126149663504
0.00031540265702
7.8852445492217E-5
1.9713222701778E-5
4.9283126335453E-6
1.2320785933717E-6
3.0801967554339E-7
7.7004920662205E-8
1.9251229943507E-8
4.8128074858766E-9
65536
131072
262144
524288
1048576
2097152
4194304
8388608
16777216
33554432
67108864
134217728
268435456
536870912
1073741824
3.14159265238659
3.14159265328899
3.14159265351459
3.14159265357099
3.14159265358509
3.14159265358862
3.1415926535895
3.14159265358972
3.14159265358977
3.14159265358979
3.14159265358979
3.14159265358979
3.14159265358979
3.14159265358979
3.14159265358979
1.2032019824915E-9
3.0080027357826E-10
7.519984634996E-11
1.8799628520583E-11
4.6993520186334E-12
1.1746159600534E-12
2.9354296771089E-13
7.327471962526E-14
1.8207657603853E-14
4.4408920985006E-15
8.8817841970013E-16
0
0
0
0
JIS-BASIC で、数値は浮動小数点 10 進数です。単精度、倍精
度の区別や、整数「型」というのはありません。10 進小数を 2 進
小数で近似するために発生する誤差がなく、0.1 を 10 回加算す
れば正確に 1 になります。たとえば、次のプログラムでは、11 回繰
り返しを実行し、最後の結果は正確に 0 になります。
(JIS-BASIC)
10 FOR x=-1 TO 0 STEP 0.1
20 PRINT x
30 NEXT x
40 END
Result
-1
-.9
-.8
-.7
-.6
-.5
-.4
-.3
-.2
-.1
0
一方、GAMBAS です。
(GAMBAS)
' Gambas module file
Public Sub Main()
Dim x As Float
For x = -1 To 0 Step 0.1
Print x
Next
Stop
End
(Result)
-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
-1.3877787807814E-16
ぐぬぬ、これでもダメか？
(GAMBAS)
' Gambas module file
Public Sub Main()
Dim x As Float
x = -1.0
x = x + 1.0E-01
x = x + 1.0E-01
x = x + 1.0E-01
x = x + 1.0E-01
x = x + 1.0E-01
x = x + 1.0E-01
x = x + 1.0E-01
x = x + 1.0E-01
x = x + 1.0E-01
x = x + 1.0E-01
Print x
Stop
End
(Resutl)
-1.3877787807814E-16
ねぇねぇ　ねぇ 今どんな気持ち?　今どんな気持ち? ∩＿＿∩　∩＿＿∩ ♪　/ ⌒　⌒丶　/⌒　⌒ 丶 ♪ /(●)　(●) ﾊｯ　(●)　(●)| ｜　(_●_) |　ﾊｯ　ﾊｯ | (_●_)　| ＿＿彡　|∪|　ミ　ﾊｯ　彡　|∪|　ミ＿＿ ＼＿＿　丶ノ>　>　<　< 丶ノ　＿＿／ /　／(_／ ♪　♪　＼_)＼　丶 /　／　＼　丶 （_⌒丶′　(フ￣つ Ｙ⌒＿) `| /ヽ|　/　丶. ｜/丶| （ﾉ　Ｕ　:｜　ー |: (ﾉ　Ｕ : >､　U(_● :/　＿ /ミノ. ｜| 三｜ |ﾟ｡ ｜| 三｜ |: (_( ￣(＿): 結果が同じ数値だけに、同程度の誤差が発生してるんだろう
なあ。(*1)
演算結果の正確さ
JIS-BASIC には数値演算の実行結果の正確さに関する規定
があります。この規定は四則演算ばかりでなく、べき乗演算や
SQR 関数などの組込み関数にも適用されます。その規定のおか
げで、べき乗演算や SQR 関数の予期せぬ誤差を心配する必要
がなくなります。
次のプログラムは 100 以下のピタゴラス数をすべて求めるも
のですが、JIS-BASIC では正しいプログラムです。しかし、規格外
の BASIC では正しい結果が得られることが保証できません。実
際、昔の N88-BASIC では正しくない結果が得られます。
(JIS-BASIC)
10 FOR x=1 TO 100
20 FOR y=x TO 100
30 LET z=SQR(x^2+y^2)
40 IF INT(z)=z THEN PRINT x,y,z
50 NEXT y
60 NEXT x
70 END
(GAMBAS)
' Gambas module file
Public Sub Main()
Dim x As Integer
Dim y As Integer
Dim z As Float
For x = 1 To 100
For y = x To 100
z = Sqr(x ^ 2 + y ^ 2)
If Int(z) = z Then
Print x, y, z
Endif
Next
Next
Stop
End
(Result)
3
4
5
12
6
8
7
24
8
15
9
12
9
40
10
24
11
60
12
16
12
35
13
84
14
48
15
20
15
36
16
30
16
63
18
24
18
80
20
21
20
48
20
99
21
28
21
72
24
32
24
45
24
70
5
13
10
25
17
15
41
26
61
20
37
85
50
25
39
34
65
30
82
29
52
101
35
75
40
51
74
25
27
28
28
30
30
32
33
33
35
36
36
39
39
40
40
40
42
45
48
48
48
51
54
56
57
60
60
60
63
65
66
69
72
75
80
60
36
45
96
40
72
60
44
56
84
48
77
52
80
42
75
96
56
60
55
64
90
68
72
90
76
63
80
91
84
72
88
92
96
100
84
65
45
53
100
50
78
68
55
65
91
60
85
65
89
58
85
104
70
75
73
80
102
85
90
106
95
87
100
109
105
97
110
115
120
125
116
うん、よし。おｋ
参考文献
・Microsoft BASIC 習熟者のための Full BASIC 入門
・白石 和夫、(仮称)十進 BASIC による JIS Full BASIC 入門の手引き、1999.1.18
・森本光生、長岡亮介、平尾淳一、奥村晴彦、数学とコンピュータ、放送大学教材、
2003.3.20
(*1)これはコンピュータの「間違い」ではなく、浮動小数点から必然的に由来する結
果というべきである。この「間違い」は決してコンピュータの欠点ではない。十進
BASIC は数学教育的観点からこのような誤差を避けるための特別処置を行なって
いる。
素数の計算
２ ch の十進 BASIC スレに先人が載せていたプログラムを
GAMBAS に移植してみました。
コンソールに数値を出力するだけですので、コマンドラインアプ
リケーションを選んでください。
最初の while ループのＣの上限値を変えてやればそこの値ま
で延々と素数をダンプします。
マシンパワーに酔い痴れてください、とのことです。
' Gambas module file
' sosu2jis.bas => sosu140916gam.bas
Public Sub Main()
Dim C As Integer
Dim K As Integer
C = 11
While C < 1000
K=3
While (Int(C / K) <> (C / K)) And (C > K * K)
K=K+2
Wend
If Int(C / K) <> (C / K) Then
Print C
End If
C=C+2
Wend
Stop
End
十進 BASIC のソースファイルも書いておきます。REM １３０行
目のＣの値を変えてやればそこの値まで延々と素数をダンプしま
す。
（JIS-BASIC）
10 OPTION ARITHMETIC RATIONAL
20 LET C=11
30 LET K=3
40 IF INT(C/K)<>(C/K) THEN
50 IF C>K*K THEN
60
LET K=K+2
70
GOTO 40
80 END IF
90 END IF
100 REM
110 IF INT(C/K)<>(C/K) THEN PRINT C;
120 LET C=C+2
130 IF C>10000 THEN STOP
1221 GOTO 30
END
二次方程式の解
まがいなりにも数値計算についてプログラミングの講釈をしよ
うというのなら、まずここからはじめなければなりませんでした。今
回は二次方程式の解です。
方法は中学校で習う、二次方程式の解そのものです。二次方
程式、Y=A×X2+B×X+C があります。定数 A,B,C を与えて、Y=0
を満たす解 X1,X2 を求めます。
ずいぶん「基本だぜ！」みたいなことを書きましたが、実は先人
が作ったソースを GAMBAS への焼き直しました。というか工夫
のしようがない。むしろこの方のウェブサイトを紹介するためにコ
ードを引用したということにしてください。他のページをご覧にな
ってみてください。この方は含蓄深いことを書いてらっしゃいます。
・ソース＆参考
教育利用を目的とした簡易な Basic インタプリタ PlainBasic プログラミング例題、
2006.08.20 版
(GAMBAS)
' Gambas module file
' PlaineBASIC ex21equation2.bas => pbasex21gam.bas
' 二次方程式を AX^2+BX+C=0 解く
Public Sub Main()
Dim A As Float
Dim B As Float
Dim C As Float
Dim D As Float
Dim X1 As Float
Dim X2 As Float
A=2
B=3
C = -3
D=B*B-4*A*C
If D < 0 Then
Print "根がありません"
Else If D = 0 Then
Print "等根です X1=X2="; - B / (A + A)
Else
D = Sqr(D)
X1 = (- B + D) / (A + A)
X2 = (- B - D) / (A + A)
Print "X1= "; X1, "X2= "; X2
Endif
Stop
End
(Result)
ex1
A=2, B=3, C=-3 の場合
X1= 0.68614066163451
X2= -2.18614066163451
ex2
A=1, B=2, C=1 の場合
等根です X1=X2=-1
ex3
A=1, B=2, C=5 の場合
根がありません
三次方程式の解
前回の続きです。これもやはり参考コードを元に GAMBAS に
焼き直しました。
三次方程式、Y=X3+A×X2+B×X+C=0 を満たす X の解を計算
します。
ここではカルダノ(Cardano,1501-1576)の公式を使います。
具体的には、3 次の対称マトリックスの固有値を計算するとき
に使うことができます。
実用例は、鉄筋コンクリート梁の設計計算のときに現われる・・
・そうです。この参考にしたウェブページを作った方は硬派なんで
すなぁ（技術者的に）。
3 つ実数解があるとき、とくに大小順に並べ直しはしてません。
・ソース＆参考
教育利用を目的とした簡易な Basic インタプリタ PlainBasic プログラミング例題、
2006.08.20 版
(GAMBAS)
' Gambas module file
' PlaineBASIC ex22equation3.txt => pbasex22gam.bas
' 三次方程式の解(カルダノの解法)
' Cardano の公式を応用して
' 三次方程式を X^3+AX^2+BX+C=0 解く
Public Sub main()
Dim P As Float
Dim A As Float
Dim B As Float
Dim C As Float
Dim D As Float
Dim Q As Float
Dim D2 As Float
Dim U1 As Float
Dim U3 As Float
Dim V1 As Float
Dim X1 As Float
Dim X2 As Float
Dim X3 As Float
Dim ANG3 As Float
Dim CC As Float
Dim SS As Float
Dim R As Float
A = -6
B = 11
C = -6
P=A*A/9-B/3
Q = -2 * A * A * A / 27 + A * B / 3 - C
D2 = Q * Q - 4 * P * P * P
D = Sqr(Abs(D2))
If D2 > 0 Then
U3 = (Q + D) / 2
U1 = (Abs(U3)) ^ (1 / 3)
If U3 < 0 Then
U1 = - U1
V1 = P / U1
X1 = U1 + V1 - A / 3
Print "実数解は一つです X1="; X1
Endif
Else If D2 = 0 Then
X1 = (4 * Q) ^ (1 / 3)
Print "等根です X1=X2=X3="; X1 - A / 3
Else
ANG3 = Atn2(D, Q) / 3
CC = Cos(ANG3)
SS = Sin(ANG3)
R = Sqr(P)
Print "三つの実根があります"
X1 = R * (CC + CC)
Print "X1= "; X1 - A / 3
X2 = R * (- CC - Sqr(3) * SS)
Print "X2= "; X2 - A / 3
X3 = R * (CC - Sqr(3) * SS)
Print "X3= "; X3 - A / 3
Endif
Stop
End
(Result)
ex1
A=-6, B=11, C=-6 の場合
三つの実根があります
X1= 3
X2= 1
X3= 2
ex2
A=-5, B=8, C=-4 の場合
実数解は一つです X1=1
ex3
A=-3, B=3, C=-1 の場合
等根です X1=X2=X3=1
方程式の解について、放送大学のテキストにわかりやすく説明
されてましたので該当ページを抜粋させていただきます。
いずれの抜粋も
・森本光生、長岡亮介、平尾淳一、奥村晴彦、数学とコンピュータ、放送大学教材、
2003.3.20
からです。
「方程式を満たす未知数 x の値(根 root, 解 solution)を求めることを方程式を解く
(solve)という。K が複素数のときは、 n 次方程式には、（重複も考慮して）n 個の根が
存在する（「代数学の基本定理」）。
1 次方程式 ax+b=0 の解は x =− b である。
a
−b±√ b2−4ac
2 次方程式 ax2+bx+c=0 の解は x =
である。
2a
2 次方程式の解から分るように、数体系が実数のみであるなら方程式の根が存
在するとは限らないが、複素数まで考慮するなら重複度を含めて 2 個の根が存在す
ることになる。
3 次方程式、4 次方程式についてはそれぞれ Cardano の解法,Ferrari の解法など
と呼ばれる代数的解法が知られている。 」
（この箇所までは部分的に省略してます）
さて、ここからが白眉です。
「Cardano の解法とは、3 次方程式
ax 3+ bx 2 + cx + d =0
を
x 3+ 3px+ q =0
という形に変形し、
x= u + v
とおくと、
( u + v )3 +3p( u + v )+ q=0
u 3+ v 3+ 3( uv + p)( u + v )+ q=0
となることから,
3
3
uv =−p かつ u + v =−q
を満たすように u, v を決めれば良い、とするものである。そのためには、
u 3 v 3 =−p3 かつ u 3+ v 3=−q
を満たす u3, v3 が、2 次方程式
t 2 + qt − p3=0
2
3
の 2 根、すなわち、 −q ±√ q + 4p であることから、
2
uv =−p
を考慮して u, v を見い出せば良いということになる。
2
3
p,q が q2+4p3>0 を満たす実数のときは、 −q ±√q + 4p が実数であるから、
2
α=
√
3
√
2
3
3 − q − √ q +4p
−q + √ q 2 +4p3
,β =
2
2
さらに、1 の原始立方根(虚数 3 乗根)の一つを ω として
{u , v }={α, β }, {αω, βω2} , {αω2 ,βω }
これより、
x =α + β ,αω + βω2 ,αω2 + βω
という 3 根を得る。
しかし,「p, q が q2+4p3>0 を満たす実数のとき」という前提が崩れたときは、この
解法は必ずしも我々が理解しているような意味での解法にはなっていない。実際、3
次方程式 x3-63x-162=0 に上の解法を形式的に適用すると、
3
3
α=√ 81+ 30 √−3 ,β= √81−30 √−3
のようなものが現れるが、
3
3
x =√ 81+ 30 √−3+ √ 81−30 √−3
は、初等的な解法によって見い出される 3 次方程式 x3-63x-162=0 の解 x=-3, -6, 9
との関係は不明である。このような事態は昔から「非還元の場合」と呼ばれてきたも
のである。
このような困難は、4 次方程式に関しても存在する。3 次方程式の解法が、もっと
も基本的な 3 次方程式 x3=a と 2 次方程式に帰着されたように、4 次方程式の解法
が、もっとも基本的な 4 次方程式 x4=a と 3 次方程式に帰着されるからである。
さらに、5 次以上の方程式については代数的な一般的解法が存在しないことが
証明されている（「Abel-Galois の定理」）。
代数学の基本定理とこれらを組み合わせると、“存在は証明されているが代数的
な方法では見つけることができない”（いわば神は知っているが人間には発見でき
ない）ということになる。 」
じゃあ 5 次方程式以上の解はどうすんの、といった感じでしょう。
そこまで来ると、方程式を解くのではなくて、解を探すという方向
に進みます。
テイラー展開
級数というのは、関数を整関数で近似しようという考え方です。
代表的な無限級数の表現方法にテイラー展開というのがあり、
exp(x)、sin(x)、cos(x)の展開式
exp(x)=1+x+x2/2!+x3/3!+x4/4!+......
sin(x)=x-x3/3!+x5/5!-x7/7!+.....
cos(x)=1-x2/2!+x4/4!-x6/6!+.....
がよく使われます。
右辺の級数は、１点 x における高次微分によって決まるもので
すから、点 x 近傍の情報から決まります。いわば局所的に決定す
るものです。左辺の意味するところは、変数 x の値を任意にとれる
ことから、いわばグローバルなものと言えるでしょう。
つまりテーラー展開というのは、微分という手段を繰り返し使
って、局所的な情報を取り出し、それを使って全体を再構成しよう
という試みなわけです。この発想の飛躍はすごいもんです。
自然対数 e を展開式と GAMBAS が持ってる関数 EXP(X)と
比較するプログラムを示します。
nmax：項数の最大値
x：変数値
x=1 の場合は自然対数の底 e を与える。
次数 n=8 で、exp(1)=2.7182539....となり、小数点以下 5 桁まで
正確な値を与える。x の値にもよるけれども、低次の展開でもよい
近似ができることがわかります。
(GAMBAS)
' Gambas module file
' taylor1gam.bas
' コンピュータ物理の功罪
' exponential(True BASIC = > GAMBAS)
Public Sub Main()
Dim nmax As Integer
Dim n As Integer
Dim x As Float
Dim t As Float
Dim s As Float
Print " nmax ? ";
Input nmax
Print " x ? ";
Input x
Print "order cal.
Let n = 1
Let t = 1.0
Let s = 1.0
Print n, s, Exp(x) - s
For n = 2 To nmax Step 1
value error"
Let t = t * x / (n - 1)
Let s = s + t
Print n, s, Exp(x) - s
Next
Print " exp(x) ="; Exp(x)
Print " exp(cal)="; s
Stop
End
(Result)
nmax ? 8
x ?1
order cal.
value error
1
1
1.71828182845905
2
2
0.71828182845905
3
2.5
0.21828182845905
4
2.66666666666667
0.05161516179238
5
2.70833333333333
0.00994849512571
6
2.71666666666667
0.00161516179238
7
2.71805555555556
0.00022627290349
8
2.71825396825397
2.7860205076724E-5
exp(x) =2.71828182845904
exp(cal)=2.71825396825397
sin の近似
振動問題では x が小さい場合、sin(x)≒x と線形近似して計算
することが結構な頻度であります。
sin(x)=x-x3/3!+x5/5!-x7/7!+.....
x<1 ではよい近似になります。項数を多くとれば x>1 でもよい
近似値を返します。
(GAMBAS)
' Gambas module file
' taylor2gam.bas
' コンピュータ物理の功罪
' sin (True BASIC = > GAMBAS)
Public Sub Main()
Dim n As Integer
Dim nmax As Integer
Dim x As Float
Dim t As Float
Dim sn As Float
Print " nmax ? ";
Input nmax
Print " x ? ";
Input x
Print "order cal.
value error"
n=2
t=x
sn = x
For n = 2 To nmax Step 2
t = - t * x ^ 2 / (n * (n + 1))
sn = sn + t
Print n, sn, Sin(x) - sn
Next
Print "sin(x) ="; Sin(x)
Print "sin(cal)="; sn
Stop
End
(Result)
nmax ? 10
x ?2
order cal.
value error
2
0.66666666666667
0.24263076015901
4
0.93333333333333
-0.02403590650765
6
0.90793650793651
0.00136091888917
8
0.90934744268078
-5.0015855094343E-5
10
0.9092961359628
1.2908628790331E-6
sin(x) =0.90929742682568
sin(cal)=0.9092961359628
GAMBAS の sin 関数(実線)と n 次のテイラー展開(x)を比べた
図を表示します。x を 0.0 から始めて For〜Next で 0.01 ずつ増加
させてます。n=13 くらまでとると x=6.28 までよく再現します。図に
してみると級数でよく近似できるもんだなぁ、と関心します。
(GAMBAS)
' Gambas class file
' taylor3gam.bas
' Taylor expansion of sin (True BASIC => GAMBAS)
Public Sub _new()
End
Public Sub Form_Open()
DrawingArea1_Draw()
End
Public Sub DrawingArea1_Draw()
Dim x As Float
Dim dx As Float
Dim t As Float
Dim sn As Float
Dim n As Integer
Dim nmax As Integer
Dim x_interval As Float
Dim y_interval As Float
Dim x_origin As Integer
Dim y_origin As Integer
Dim x1 As Integer
Dim y1 As Integer
Dim x2 As Integer
Dim y2 As Integer
Dim xp As Integer
Dim yp As Integer
'
'
'
'
'
Print " nmax(<15) ? "
Input nmax
x = 0.0
dx = 0.01
Initialize the DrawingArea object
DrawingArea1.Clear
データを描画
draw.Begin(DrawingArea1)
縦軸と横軸の目盛りを決める。表示する項目の数、データの最大値によって変わる。
x_interval = DrawingArea1.Width / (6.28 - 0.0)
y_interval = DrawingArea1.Height / (1.0 - (-1.0))
Define the origin for the chart
x_origin = Int(Abs(0.0) * x_interval)
y_origin = DrawingArea1.Height - Int(Abs(-1.0) * y_interval)
x1 = x_origin
y1 = y_origin
For x = 0 To 6.28 Step dx
t=x
sn = x
x2 = x_origin + Int(x * x_interval)
y2 = y_origin - Int(Sin(x) * y_interval)
draw.line(x1, y1, x2, y2)
For n = 2 To nmax + 1 Step 2
xp = x_origin + Int(x * x_interval)
yp = y_origin - Int(sn * y_interval)
draw.point(xp, yp)
t = - t * x ^ 2 / (n * (n + 1))
sn = sn + t
Next
x1 = x2
y1 = y2
Next
Finish with the DrawingArea object
draw.End
End
(Result)
参考文献
・十進 BASIC ではじめる BASIC、講談社ブルーバックス
・和田正信著、有馬朗人、大槻義彦編、コンピュータ物理の功罪、物理数学
OnePoint10、共立出版
で、この参考文献ですが。『コンピュータ物理の功罪』と言いま
すと、やってることを否定すんのかよッ！と煽ってる題名のように感
じますが、中身はちゃんとしたことが書かれてます。
このシリーズの編者のお二人はかなり有名な方々なので特に
説明する必要はないと思います。この本が出版された当時は、32
ビット PC が普及しはじめた時代でして、大学の研究室でも PC で
数値計算させるようになってきた頃だと思います。そういうパーソ
ナルコンピュータが普及するのを見て、やはり物理屋であるお二
人は、「数値計算だけまとめてわかったふりしてんじゃないぞ」と
釘を差すつもりもあったんじゃないでしょうか。ただの妄想ですが。
まあ世間でよく聞く慣用句に、「そんなこと言っても、理屈と現
実は違うんだ！」と言うセリフがありますが、これはまさにそういう
ことです。理屈と現実が違うように見えてるのは、そう言う人。その
方の世界の見方を言い表してるんだなと思うわけです。そんなこ
たない。理屈と現実が違うのなら、なぜ違うのかを考えるのが人
間というものです。
ニュートン法
ここまでくれば、だいぶ数値計算してるっぽくなりますよね。
例によって他人様が作ったコードの GAMBAS への焼き直し
です。根性なしです。「根なし」だけに。
一般的な問題の応用には、解きたい関数本体 Y、その一回微
分した式 dY を差し換えればよいです。
Newton 法では、初期値から出発して補正を繰り返して一つの
解を得ます。二つ以上の解があるときは、あらかじめ、どのあたり
に解があるか初期値を入れ直さなければなりません。
また、繰り返し計算を止める条件として、精度を決めておかな
ければなりません。もし極大または極小値があると無限ループに
入る危険があります。これを避けるため、For-Next 文で 50 回の繰
り返しで計算を停止するようにしてあります。
・ソース＆参考
教育利用を目的とした簡易な Basic インタプリタ PlainBasic プログラミング例題、
2006.08.20 版
(GAMBAS)
' Gambas module file
' 三次方程式の解(ニュートン法)
' ex23Newton.bas => pbasex23gam.bas
' Nweton の逐次近似法を使って
' 三次方程式 X^3+AX^2+BX+C=0 を解く
Public Sub Main()
Dim A As Float
Dim B As Float
Dim C As Float
Dim XN As Float
Dim ERR As Float
Dim X As Float
Dim Y As Float
Dim DY As Float
Dim Count As Integer
A = -6 ' 例題です
B = 11
C = -6
XN = 0 ' 初期値です
' XN=10 ' 別の初期値です
ERR = 1.0E-15 ' 精度設定です
X = XN ' 初期値です
Y = ((X + A) * X + B) * X + C
DY = (3 * X + 2 * A) * X + B
XN = X - Y / DY
While ((Abs(XN - X) > ERR) And (Count < 50))
X = XN
Y = ((X + A) * X + B) * X + C
DY = (3 * X + 2 * A) * X + B
XN = X - Y / DY
Count = Count + 1
Print Count, "X="; X, "XN="; XN, "XN-X="; XN - X
Wend
Stop
End
(Result)
例題
A = -6
B = 11
C = -6
初期値 XN = 0 の場合
1 X=0.54545454545455 XN=0.84895321062245 XN-X=0.30349866516791
2 X=0.84895321062245 XN=0.9746740710233 XN-X=0.12572086040085
3 X=0.9746740710233 XN=0.99909154805695 XN-X=0.02441747703365
4 X=0.99909154805695 XN=0.99999876469105 XN-X=0.00090721663411
5 X=0.99999876469105 XN=0.99999999999771 XN-X=1.2353066561888E-6
6 X=0.99999999999771 XN=1
XN-X=2.2888357875672E-12
7 X=1
XN=1
XN-X=0
別の初期値 XN = 10 の場合
1 X=7.36125654450262 XN=5.61610701395443 XN-X=-1.74514953054819
2 X=5.61610701395443 XN=4.47379897271395 XN-X=-1.14230804124047
3 X=4.47379897271395 XN=3.74420450714993 XN-X=-0.72959446556402
4 X=3.74420450714993 XN=3.30588643643379 XN-X=-0.43831807071614
5 X=3.30588643643379 XN=3.08210386279587 XN-X=-0.22378257363792
6 X=3.08210386279587 XN=3.00848840787698 XN-X=-0.07361545491889
7 X=3.00848840787698 XN=3.00010598093566 XN-X=-0.00838242694132
8 X=3.00010598093566 XN=3.00000001684377 XN-X=-0.00010596409189
9 X=3.00000001684377 XN=3
XN-X=-1.684377126665E-8
10 X=3
XN=3
XN-X=-8.8817841970013E-16
二元連立一次方程式
行列の計算です。逆行列の計算に至って、パソコンやスマフォ
を計算機として使ってる感が出てくると思います。
これは典型的な二元連立一次方程式を解く問題で、二つの直
線の交点を求めます。やっぱりこれも中学校でならう範囲ですよ
ね。この例では、6 個の定数を入力データとします。リストの変数
値を他の数値に差しかえると、別の計算ができます。
まあ GAMBAS に限らず、プログラミング言語のほとんどが行
列記号そのものを使うという文法がないので、代数式の形で扱う
ことになります。言語（コンパイラ）という括りではなく、数式処理
アプリという分類では FreeMat とか Octave というフリーウェアが
世に出ていますので、そちらもご参考にしてください。
・ソース＆参考
教育利用を目的とした簡易な Basic インタプリタ PlainBasic プログラミング例
題、2006.08.20 版
(GAMBAS)
' Gambas module file
' PlainBASIC ex24lineareq2.bas => pbasex24gam.bas
' 二元連立一次方程式を解く
' 2 直線の交点を求める問題は二元一次連立方程式を解きます
' A*X+B*Y=P
' C*X+D*Y=Q
Public Sub Main()
Dim A As Float
Dim B As Float
Dim P As Float
Dim C As Float
Dim D As Float
Dim Q As Float
Dim DET As Float
Dim X As Float
Dim Y As Float
A=1
B=2
P=5
C=3
D = -1
Q=1
DET = A * D - C * B
If Abs(DET) < 1.0E-15 Then
Print "2 直線は平行です"
Else
X = (P * D - Q * B) / DET
Y = (Q * A - P * C) / DET
Print "X= "; X
Print "Y= "; Y
Endif
Stop
End
(ex1)
A = 1, B = 2, P = 5
C = 3, D = -1, Q = 1
(Reslut1)
X= 1
Y= 2
例：平行線の場合
DET がゼロになると、ゼロで割ることになるので DET の値が GAMBAS の有効
桁数で 1.0E-15 より小さければ平行線と判定する。
(ex2)
A = 1, B = 2, P = 5
C = 1, D = 2, Q = 1
(Reslut2)
DET = 0
2 直線は平行です
(ex3)
A = 1.0 + 1.0E-15, B = 2, P = 5
C = 1, D = 2, Q = 1
(Result3)
DET = 2.2204460492503E-15
X= 3.6028797018964E+15
Y= -1.8014398509482E+15
(ex4)
A = 1.0 + 1.0E-16, B = 2, P = 5
C = 1, D = 2, Q = 1
(Result4)
DET = 0
2 直線は平行です
逆行列の計算
逆行列を求めることができれば、連立方程式も解けるよね。と
いうわけで、数値計算で逆行列を求める方法は、ほんと良く使い
ます。
逆行列の求め方となると、はじめはクラメルの方法を教わるこ
とになります。けれど、N 次の逆行列の求め方となると収拾がつ
かなくなるでしょう。
大学教養過程では『基礎線形代数』とか『線形代数の基礎』
という講義でマシな方法を習います。早い話がいかに係数を減ら
せるかにかかってます。実際に教科書の演習問題を解いてみると
記号の操作が良く理解できると思います。ここでは数値計算でよ
く使われるガウスの方法を紹介します。
（１）ガウス・ジョルダン法
ガウス・ジョルダン法は「吐き出し法」ともいわれ、次元が大き
い場合、クラメル法に比べて演算回数が少なくてすむため数値
解法としてよく使われます。N 元の場合、乗除算では
N(N2-3N-1)/3 回で、N=5 で 65 回、N=10 でも 430 回という感じで
す。
まとめると、N 元連立１次方程式の場合、以下のことを実行しま
す。
1）ai,k←ai,k/ai,i を i=1,2,...,N について繰り返し、k=i+1,...,N-1 まで
繰り返す。
2）aj,N+1←aj,N+1-aj,i・ai,k を j=1,...,N(j≠i)まで繰り返し、k=i+1,...,N-1 ま
で繰り返す。
3）xj=aj,N+1 で解を求める。
『コンピュータ物理の功罪』(p.41)より
(例題)次のような 4 元連立方程式を解きます。
3 −4 1 6 x 0
22
−2 10 −5 3 x 1 = 15
1
2 −2 8 x 2
31
1 −5 7 7 x 3
40
(
)( ) ( )
(GAMBAS)
' Gambas module file
Public Sub Main()
' GaussJordanGam.bas
' TrueBASIC => GAMBAS
' Linear Equation(Gauss-Jordan method)
Dim i As Integer
Dim j As Integer
Dim k As Integer
Dim n As Integer = 4
Dim a As New Float[n, n+1]
a[0, 0] = 3.0
a[0, 1] = -4.0
a[0, 2] = 1.0
a[0, 3] = 6.0
a[0, 4] = 22.0
a[1, 0] = -2.0
a[1, 1] = 10.0
a[1, 2] = -5.0
a[1, 3] = 3.0
a[1, 4] = 15.0
a[2, 0] = 1.0
a[2, 1] = 2.0
a[2, 2] = -2.0
a[2, 3] = 8.0
a[2, 4] = 31.0
a[3, 0] = 1.0
a[3, 1] = -5.0
a[3, 2] = 7.0
a[3, 3] = 7.0
a[3, 4] = 40.0
k=0
While k < n
For j = k + 1 To n
Let a[k, j] = a[k, j] / a[k, k]
Next
For i = 0 To n - 1
If i <> k Then
For j = k + 1 To n
a[i, j] = a[i, j] - a[i, k] * a[k, j]
Next
Endif
Next
k=k+1
Wend
For i = 0 To n - 1
Print " X"; i; " = "; a[i, n]
Next
Stop
End
(Result)
X0 = 1
X1 = 2
X2 = 3
X3 = 4
（２）ガウスの消去法
未知数を１個ずつ消去し、最終的に 0 でない係数が 1 個だけ
の方程式を作り（前進消去法）、その解を前の式に順次代入して
（後進代入）すべての未知数を求める方法です。ガウスの消去法
はガウス・ジョルダン法より計算手順が複雑なんですが、必要な
演算回数が 6 割くらい少なくて数値計算向きの方法です。
1）ak,j←ak,j/ak,k を k=1,...,N-1、j=k+1,...,N まで繰り返す。
2）ai,j←ai,j-ai,k・ak,j、ai,N+1←ai,N+1-ai,N+1・ak,k を i=k+1,...,N、j=k,...,N ま
で繰り返す。
3）xN←aN,N+1/aN,N
4）xk←ak,N+1-ak,jxj を k=N-1,...,1 と減らしながら解を求める。
ガウス・ジョルダン法の例題と同じ行列式を解きます。結果は
同じになるはずですね。
(GAMBAS)
' Gambas module file
Public Sub Main()
' GaussmethodGam.bas
' TrueBASIC => GAMBAS
' Linear Equation-Gauss method
Dim i As Integer
Dim j As Integer
Dim k As Integer
Dim n As Integer = 4
Dim p As Float
Dim s As Float
Dim a As New Float[n, n + 1]
Dim x As New Float[n]
a[0, 0] = 3.0
a[0, 1] = -4.0
a[0, 2] = 1.0
a[0, 3] = 6.0
a[0, 4] = 22.0
a[1, 0] = -2.0
a[1, 1] = 10.0
a[1, 2] = -5.0
a[1, 3] = 3.0
a[1, 4] = 15.0
a[2, 0] = 1.0
a[2, 1] = 2.0
a[2, 2] = -2.0
a[2, 3] = 8.0
a[2, 4] = 31.0
a[3, 0] = 1.0
a[3, 1] = -5.0
a[3, 2] = 7.0
a[3, 3] = 7.0
a[3, 4] = 40.0
For k = 0 To n - 2
For i = k + 1 To n - 1
p = a[i, k] / a[k, k]
For j = k To n
a[i, j] = a[i, j] - p * a[k, j]
Next
Next
Next
x[n - 1] = a[n - 1, n] / a[n - 1, n - 1]
For k = n - 2 To 0 Step -1
s = a[k, n]
For j = k + 1 To n - 1
s = s - a[k, j] * x[j]
Next
x[k] = s / a[k, k]
Next
For i = 0 To n - 1
Print " x"; i; " = "; x[i]
Next
Stop
End
(Result1)
x0 = 0.99999999999999
x1 = 1.99999999999999
x2 = 2.99999999999999
x3 = 4
違う参考文献から持って来たソースも載せておきます。同じ計
算を違うソースで比較してみると、コードのポイントがわかりやす
いと思いますので。こちらの方は、元ネタが ExcelVBA マクロだっ
たのを GAMBAS に焼き直しました。 前の例題と同じ式を解きま
す。
(GAMBAS)
' Gambas module file
' gaussvbagam.bas
' エクセルプログラム => GAMBAS
' プログラム 1.1 ガウスの消去法
' 表 1.1 ガウスの消去法プログラムリスト
' 元数...2 元 1 次方程式であれば 2 となり 4 元 1 次方程式であれば 4 になる
Public n As Integer = 4
Public A As New Float[n, n]
Public B As New Float[n]
Public X As New Float[n]
Public Sub Main()
Dim i As Integer
A[0, 0] = 3.0
A[0, 1] = -4.0
A[0, 2] = 1.0
A[0, 3] = 6.0
A[1, 0] = -2.0
A[1, 1] = 10.0
A[1, 2] = -5.0
A[1, 3] = 3.0
A[2, 0] = 1.0
A[2, 1] = 2.0
A[2, 2] = -2.0
A[2, 3] = 8.0
A[3, 0] = 1.0
A[3, 1] = -5.0
A[3, 2] = 7.0
A[3, 3] = 7.0
B[0] = 22.0
B[1] = 15.0
B[2] = 31.0
B[3] = 40.0
Calculation(n)
' 解
For i = 0 To n - 1
Print " X"; i; " = "; X[i]
Next
Stop
End
' ***** ガウスの消去法
Sub Calculation(n As Integer)
Dim PVT As Float
Dim DIAG As Float
Dim AM As Float
Dim S As Float
Dim h As Integer
Dim i As Integer
Dim j As Integer
Dim k As Integer
PVT = 0.00005
' ***** 前進消去
For h = 0 To n - 1
If (Abs(A[h, h]) < PVT) Then
Print "ピボットが" & h & "番目の消去で小さくなり過ぎました"
Stop
Endif
DIAG = 1.0 / A[h, h]
For j = h + 1 To n - 1
AM = A[j, h] * DIAG
For k = h + 1 To n - 1
A[j, k] = A[j, k] - AM * A[h, k]
Next
B[j] = B[j] - AM * B[h]
Next
For k = h + 1 To n - 1
A[h, k] = A[h, k] * DIAG
Next
B[h] = B[h] * DIAG
Next
' ***** 後退代入
X[n - 1] = B[n - 1] / A[n - 1, n - 1]
For j = n - 1 To 0 Step -1
S = 0.0
For k = j + 1 To n - 1
S = S + A[j, k] * X[k]
Next
X[j] = B[j] - S
Next
End
(Result2)
X0 = 0.99999999999999
X1 = 2
X2 = 3
X3 = 4
（３）コレスキー法
おそらく、ポケコンジャーナルに掲載されていたプログラムを
GAMBAS に書き直してみました。手元に残っていたのがコピー
の切れ端のため初出は不確かです。
今回のプログラムのポイントは「逆行列を求める」というテーマ
を逸脱したところにありまして、プログラミングそのもののクール
さに注目して頂きたいところです。逆行列を計算するという連載
は、このコードを読んでいただきたいがため続けたと言っても過
言ではないでしょう。
ポケコンという機能が制限された環境で実に工夫してプログラ
ムしています。この頃のポケコンの BASIC では、変数名が A から
Z までのアルファベット 26 文字。そこから配列を使うとなると、そ
の文変数の数が制限されるというタイトなものだったんです。しか
も表示は１行だけ。兎にも角にも、それはメモリ容量に起因した制
約なんです。このプログラムのソースを読むと、キレーイにアルフ
ァベット 26 文字を使い切ってる！
GAMBAS に翻訳した都合上、もともと 3 つのプログラムに分
割されていたソースを１つにまとめています。GAMBAS を使えば
配列の次元を複数とるのはわけのない話ですが、このプログラム
が作られた当時の時代背景を想像して、そのダンディズムを味わ
っていただきたい。
(例題)次のような 3 元連立方程式を解きます。
(
1 2
3 x 1 −74
2 20 26 x 2 = 52
3 26 70 x 3
100
)( ) ( )
(GAMBAS)
' Gambas module file
' cholesky0216gam.bas
' 2014/02/15
' program title : コレスキー法
' PJ かな？あなたの仕事をスピーディにお手伝い（データ処理用）
' coding : 前田修三郎 氏
Public DAT As New Float[26]
Public Sub Main()
Dim A As Integer
Print "- 元の数を入力してください"
Print " ? = ";
A=3
Input A
AA1(A)
AA2(A)
AA3(A)
Stop
End
Sub AA1(A As Integer)
Dim B As Integer
Dim D As Integer
Dim G As Integer
Dim H As Integer
Dim K As Integer
K=A
B=A-1
Print "- 左辺 要素 a(G 行,H 列)の数値入力"
For G = 1 To A
D=A-B
For H = D To A
Print "("; G; ","; H; ") = ";
Input DAT[K]
K=K+1
Next
B=B-1
Next
End
Sub AA2(A As Integer)
Dim B As Integer
Dim D As Integer
Dim E As Float
Dim F As Integer
Dim G As Integer
Dim H As Integer
Dim I As Integer
Dim J As Integer
Dim K As Integer
K=A
DAT[K] = Sqr(DAT[K])
E = DAT[K]
K=K+1
For G = 1 To A - 1
DAT[K] = DAT[K] / E
K=K+1
Next
B=A-1
For G = 0 To A - 3
GOSUB500(K, B, G)
E = DAT[K]
K=K+1
For H = 1 To B - 1
F=K-B
J=0
For I = 1 To A - B
J = J + DAT[F] * DAT[F - H]
F=F-B-I
Next
DAT[K] = (DAT[K] - J) / E
K=K+1
Next
B=B-1
Next
G=A-2
GOSUB500(K, B, G)
' コレスキー分解の結果
K=A
B=A-1
Print "- コレスキー分解の結果"
For G = 1 To A
D=A-B
For H = D To A
Print "("; G; ","; H; ")="; DAT[K]
K=K+1
Next
B=B-1
Next
End
Sub GOSUB500(K As Integer, B As Integer, G As Integer)
Dim E As Integer
Dim F As Integer
Dim H As Integer
F=K-B
E=0
For H = 1 To G + 1
E = E + DAT[F] * DAT[F]
F=F-B-H
Next
DAT[K] = Sqr(DAT[K] - E)
End
Sub AA3(A As Integer)
Dim B As Integer
Dim C As Integer
Dim F As Integer
Dim G As Integer
Dim H As Integer
Dim J As Integer
Dim K As Integer
Dim L As Float
Print "- 右辺入力"
For G = 1 To A
Print "("; G; ") = ";
C=G-1
Input DAT[C]
Next
B=A-1
L = L / DAT[A]
For G = 1 To A - 1
C=0
K=C+A+G
J=0
For H = 1 To A - B
J = J + DAT[K] * DAT[C]
K=K+A-H
C=C+1
Next
B=B-1
K = A + G * A - (G * G - G) / 2
C=G
DAT[C] = (DAT[C] - J) / DAT[K]
Next
C=A-1
K = C + (A * A + A) / 2
DAT[C] = DAT[C] / DAT[K]
B=A-1
For G = 2 To A
C=A-1
J=0
F=K-1
For H = 1 To A - B
J = J + DAT[F] * DAT[C]
C=C-1
F=F-1
Next
K=K-G
C=A-G
DAT[C] = (DAT[C] - J) / DAT[K]
B=B-1
Next
Print "-Result"
For G = 1 To A
C=G-1
Print G; " = "; DAT[C]
Next
End
(Result)
- 元の数を入力してください
?=3
- 左辺 要素 a(G 行,H 列)の数値入力
(1,1) = 1
(1,2) = 2
(1,3) = 3
(2,2) = 20
(2,3) = 26
(3,3) = 70
- コレスキー分解の結果
(1,1)=1
(1,2)=2
(1,3)=3
(2,2)=4
(2,3)=5
(3,3)=6
- 右辺入力
(1) = -74
(2) = 52
(3) = 100
-Result
1 = -100
2 = 10
3=2
〜〜掲載オリジナルの解説〜〜
あなたの仕事をスピーディにお手伝い（データ処理用）
coding : 前田修三郎 氏
行列式でお悩みの方、ちょっと聞いてください（ガウスの消去法の改良版）
複雑な行列式を前に、頭を抱えているあなた。でも、よーく見てください。あなたが
計算しようとしている連立方程式の係数行列は、対称かつ正定値ではありませんか
？もしそうなら、とてもラッキー。このプログラムがあなたに代わって、コレスキー法を
用いて解を求めます。また、コレスキー法の持つ特性（掃き出し法と比べて演算回数
が少なく、対称性を利用して三角行列を解く）から、ここでは１５元までの係数行列を
計算することができます。
なお、このプログラムでは、１つのプログラムを”AA1"”AA2"”AA3"の３つに分割し
ています。CHAIN 文によって続けることになっているので、あらかじめ、３つのプログ
ラムをカセットテープにセーブしておくことが必要です。
使い方
たとえば、次のような３元連立方程式がある場合、
（１）（２）＝（３）
ただし
３つのプログラムの内容は、次のとおりです。
プログラム　内容
AA1　左辺（１）右上半分の三角行列の入力および確認・訂正
AA2　（１）のコレスキー分解・結果出力
AA3　右辺（３）の入力および演算・結果（２）出力
(GAMBAS では、AA1,AA2,AA3 を Sub プロシージャとして引数を受け渡しながら
順番に処理して行ってます。)
AA1 プログラム
１.AA1 のプログラムを開始します。
２.最初に方程式の元数を入力します。次に、表示にしたがって、左辺（１）の値を入力
していきます。表示は（例）を参照してください。
（例）”１　１”→行列の１行１列の要素 a11
３．左辺の値の入力が終わると、AA2 で、入力した値の確認を行います。”（１、１）＝１
”というふうに表示されるので、その値に誤りがなければ、[Enter]を押し、”OK THEN
INPUT 1"に１を入力してください。誤りがあった場合は、[Enter]を押して”OK THEN
INPUT 1"のところで１以外の数を入力します。続いて、"TEISEI="と聞いてくるので、
正しい値を入力してください。（以下、くり返します）
AA2 プログラム
4.確認が終わると、”AA2”のプログラムがセーブされたテープを読みこめる状態に
セットし、プログラムをロードします。
５.ロードが完了すると、自動的に計算を始めます。計算が終了すると、ピー音が３回
鳴ります。コレスキー分解の結果が表示されます。
AA3 プログラム
６.”AA3”のプログラムがセーブされたデータを読みこめる状態にセットします。
７.ロードが完了したら、右辺の値を表示にしたがい、順に入力します。
８.計算を開始し、結果（２）を表示していきます。（以下、順次入力）
９.再び解を表示させたい時は、そう入力します。
１０.さらに、左辺（１）が同一で、右辺（３）のみ異なる方程式を解きたい場合は、RUN
を入力し、７以下の操作をくり返します。
参考（計算内容など）
AA2 の左辺計算
まず C11=√ a11 .
次で i=2,3,....,n に対し Ci1=ai1 / C11 .
√
k −1
k=2,3,...,n の項に Ckk = a kk− ∑ C2kj .
j=1
k=n ならば、ここで終わり。
k −1
k<n ならば i=k+1,...,n に対し n-1 まで Cik =( aik − ∑ Cij Ckj )/ Ckk .
j=1
√
n−1
最後に Cnn = a nn− ∑ C2nj .
j=1
AA3 の右辺計算
まず y1= b1 / C11 .
n−1
k=2,3,...,n の順に yk =( bk − ∑ Cki y i )/C kk （前進消去）。
i =1
k=n-1,...,2,1 の順に x k =( y k −
となる。
数値積分
（１） 区分求積法
n
∑
j= k +1
Cik X i)/ C kk （後進代入）。
定積分は面積と関連づけるとたやすく理解できます。関数のグ
ラフと x 軸によって囲まれた面積が定積分であるというわけです。
曲線によって囲まれた図形を面積をどうやって計算するかとい
いますと、この図形を細長い長方形の束であると見なして近似す
るわけです。
積分の教科書には大抵そういう模式図が描かれてると思いま
す。では計算機で積分値 ∫f(x)dx を計算してみましょう。
x=a から x=b まで N 等分して、その幅を v=(b-a)/n とすると、 面
積 s は、
s=f(a)v+f(a+h)v+ … + f(a+(n-1)v)v
で近似される。分割数 N を増やしていけば、より正確な積分値に
近づくはずです。
サンプルプログラムは、y=x2 の 0 から 1 までの面積を求めるプ
ログラムです。関数の式は Func 関数に書いてあるとおりです。
For 文の Step 0.01 が 100 等分することを意味しています。これ
を 1000 とか 10000 に増やせば、より細かく分割されるので時間
はかかるけども厳密解に近づくはず、という寸法です。
参考文献
・木村良夫、パソコンを遊ぶ簡単プログラミング、コンピュータを自由に操る「十進
BASIC」入門、講談社ブルーバックス、2003.01.20
・大学院分子機能薬学専攻、創薬基盤分子設計学講座
情報処理教育用 BASIC 言語入門（UNIX(Linux)実習前の BASIC プログラミン
グ自習資料）
(GAMBAS)
' Gambas module file
' kubun9gam.bas
' 区分求積法
Public Sub Main()
Dim i As Float
Dim l As Float
Dim s As Float
Dim t As Float
t=0
For i = 0 To 1 Step 0.01
l = Fnc(i)
t = 0.01 * l
s=s+t
Next
Print s
Stop
End
Function Fnc(x As Float) As Float
Return x * x
End
（２）シンプソンの積分則
数値積分の公式としてよく知られている公式に次のようなのが
あります。
・ニュートン・コーツ（Newton-Coutes）系の公式
直線近似：台形則、中点則
２次曲線近似：シンプソン則
・チェビシェフ（Chebyshev）の公式
（区間を不等間隔で重みを一定にする方法）
・ガウス（Gauss）の公式
（区間を不等間隔で重みも一定でない方法。もっとも精度が高
い方法）
シンプソンの積分則のサンプルプログラムを提示いたします。x
の関数 f(x)の区間[a,b]の定積分を求めます。やってることは区分
求積法の発展バージョンです。お約束ですが、f(x)は区間[a,b]で
連続で、区間内で求めうるものとします。
f(x)を直接積分するのが難しいとき、f(x)を別の積分可能な関
数 g(x)に置き換えることで、f(x)の近似積分を求める手段を取りま
す。数値積分を使おうとするからには、厳密に求めるのが困難だ
からなわけで、それじゃ簡単な近似に置き換えようという考え方
です。そういう発想の求積法をニュートン・コーツ系の公式という
そうです。
それ系の公式で、区間[a,b]の両端点 a,b とその中間
x=(a+b)/2 で、f(x)と一致するような 2 次の多項式 g(x)を使うのを
シンプソン則といいます。たくらんでることは、なんとなくイメージ
つくと思います。
微小区間 h をさらに 2 つの区間にして積分を行うわけですの
で
h=(b-a)/2n
x(i)=a+i・h
（ここで、i=0,1,...,2n）
とおくと、x(0)=a,x(2n)=b となるように工夫する。n を大きくするほ
ど細かく分割することになるので精度がよくなります。
例えば N=10 のとき f ( x )=4 /(1+ x ) を 0 から 1 まで積分すると、
正確には 3.14159...（=PI）となります。ちなみにこれは x=tanθ と
置換積分すると解析的に求められます。
2
(GAMBAS) simpson1gam.bas
' Gambas module file
' simpsongam.bas
' Simpson Integral Method (TrueBASIC => GAMBAS)
Public Sub Main()
Dim i As Integer
Dim j As Integer
Dim m As Integer
Dim n As Integer
Dim x As Float
Dim a As Float
Dim b As Float
Dim h As Float
Dim s As Float
Dim s1 As Float
Dim s2 As Float
Print "Simpson Method"
Print " - interval [a,b] "
Print " a ? ";
Input a
Print " b ? ";
Input b
Print " - division m"
Print " m ? ";
Input m
n=2*m
h = (b - a) / n
s1 = 0.0
s2 = 0.0
s = 0.0
For i = 1 To n - 1 Step 2
x=a+i*h
s1 = s1 + fx(x)
Next
For j = 2 To n - 2 Step 2
x=a+j*h
s2 = s2 + fx(x)
Next
s = h * (fx(a) + fx(b) + 4 * s1 + 2 * s2) / 3
Print " cal.value", "exact.value", "error"
Print s, Pi, s - Pi
Stop
End
Function fx(x As Float) As Float
Dim fx As Float
fx = 4 / (1 + x ^ 2)
Return fx
End
(Result)
Simpson Method
- interval [a,b]
a?0
b?1
- division m
m ? 10
cal.value
exact.value error
3.14159265296979
3.14159265358979
また、関数 g(x)を工夫して
h=(b-a)/n
-6.200076008156E-10
x(i)=a+i・h
(ここで、i=0,...,n)
とおくと、n を奇数偶数に関係なくとれるので少し簡単になります。
(GAMBAS) simpson2gam.bas
' Gambas module file
' simpson2gam.bas
' Integral(Simpson method) (True BASIC => GAMBAS)
Public Sub Main()
Dim a As Float
Dim h As Float
Dim n As Integer
Dim b As Float
Dim s As Float
Print "Simpson method"
Print " - interval [a. b]"
Print " a ? ";
Input a
Print " b ? ";
Input b
Print " - division n"
Print " n ? ";
Input n
h = (b - a) / n
Print "cal.value", "exact.value", "error"
s = Smp3(a, h, n)
Print s, Pi, s - Pi
Stop
End
Function fx(x As Float) As Float
Dim fx As Float
fx = 4 / (1 + x * x)
Return fx
End
Function Smp3(a As Float, h As Float, n As Integer) As Float
Dim s As Float
Dim x As Float
Dim i As Integer
s=0
For i = 0 To n - 1
x=a+i*h
s = s + fx(x) + 4 * fx(x + h / 2) + fx(x + h)
Next
s=h*s/6
Return s
End
(Result)
Simpson method
- interval [a. b]
a?0
b?1
- division n
n ? 10
cal.value
exact.value error
3.14159265296979
3.14159265358979
-6.2000715672639E-10
まあ私にとってのシンプソンといえば『シンプソンズ』（黄色い
人類のでる漫画）なんですが、外国人の名前の付いたなんか憎
ったらしい公式なのかなぁ、と思いきや。わりと素直な発想で好
感がもてますね。
参考文献
・和田正信著、有馬朗人、大槻義彦編、コンピュータ物理の功罪、物理数学
OnePoint10、共立出版
・河西朝雄、VisualBasic によるはじめてのアルゴリズム入門、技術評論社、平成 11
年 6 月 10 日
（３）ガウスの積分則
区分求積法ではあたりまえのように考えてた、積分点が等間
隔であるというお約束はおしまいです。
あと積分区間の正規化を行います。今まで積分区間を a から
b までとしてましたが、それを-1 から+1 までのように規格化しま
す。そうすると問題を一般化できて楽なんですね。
また一般に、任意の関数 f(x) は n 次のルジャンドル多項式と
N-1 次以下の関数 Q,R で表される。
f ( x)=Q( x)⋅P N + R(x )
ここで Pi(x) は区間[-1,1]において次の式を満たす直交したル
ジャンドル多項式である。
1
∫ P i ( x)⋅x k dx=0 (k=0,1,2,・・・,N-1)
−1
Pi (x)=0 (i=1,2,・・・,N-1)
上の 2 つの条件を満たす i 個の解が区間[-1,1]の分点である。
この分点はルジャンドル多項式のゼロ点である。
N-1 次またはそれより低い次数の関数 R が正確に積分される
ように wi を連立した線形方程式を満たすように選ぶ。ルジャンド
ル多項式の性質を利用して
w i=
2(1−xi )2
[N P N−1 ( x i)] 2
を得る。したがってガウスの積分公式は
∫ f ( x)dx= 12 (b−a)⋅∑ w f (x )
b
N
a
i=1
i
i
で与えられる。
サンプルプログラムは、N=3 or 5 の場合の精度を
f ( x )=4 /(1+ x ) を用い、積分区間 0 から 1 まで積分して調べます。
2
(GAMBAS)
' Gambas module file
' gaussintegam.bas
' Gauss Integral (True BASIC => GAMBAS)
Public Sub Main()
Dim t As New Float[6]
Dim c As New Float[6]
Dim dx As Float
Dim hx As Float
Dim ss As Float
Dim s As Float
Dim x0 As Float
Dim x As Float
Dim a As Float
Dim b As Float
Dim i As Integer
Dim k As Integer
Dim m As Integer
Dim n As Integer
Print " Nth Gauss Integral(N=3 or 5)"
Print " interval [a, b] "
Print " a ? ";
Input a
Print " b ? ";
Input b
Print " division m ? ";
Input m
Print " degree N(3 or 5) ? ";
Input n
dx = (b - a) / m
hx = dx / 2.0
ss = 0.0
t[1] = -0.7745966692
t[2] = 0.0
t[3] = 0.7745966692
t[4] = 0.0
t[5] = 0.0
c[1] = 0.5555555556
c[2] = 0.8888888889
c[3] = 0.5555555556
c[4] = 0.0
c[5] = 0.0
If n <> 3 Then
t[1] = -0.9061798459
t[2] = -0.5384693101
t[3] = 0.0
t[4] = 0.5384693101
t[5] = 0.9061798459
c[1] = 0.2369268851
c[2] = 0.4786286705
c[3] = 0.5688888889
c[4] = 0.4786286705
c[5] = 0.2369268851
End If
For k = 0 To m - 1
x0 = a + k * dx + hx
For i = 1 To n
x = x0 + t[i] * hx
ss = ss + c[i] * f(x)
Next
Next
s = ss * hx
Print "cal.value", "exact.value", "error"
Print s, Pi, s - Pi
Stop
End
Function f(x As Float) As Float
Dim f As Float
f = 4.0 / (1 + x * x)
Return f
End
(Result)
Nth Gauss Integral(N=3 or 5)
interval [a, b]
a?0
b?1
division m ? 10
degree N(3 or 5) ? 3
cal.value
exact.value error
3.14159265371718
3.14159265358979
1.2738432531023E-10
分割数 m を適当に 1 からあげてみると、次数が低くてもシンプ
ソンの公式よりも良い結果が出てくる。
積分区間が大きい場合は区間を m 個に分割して、その小区間
にガウス積分を適用し全体の和をとると、次数 N を小さくしても
精度のよい積分値を得ることができる。
参考文献
・和田正信著、有馬朗人、大槻義彦編、コンピュータ物理の功罪、物理数学
OnePoint10、共立出版
常微分程式の数値解法
関数 y=f(x)と、これを微分したもの y'=f'(x)、あるいはこれをさら
に微分したもの y''=f''(x)などの間の関係を表す式を常微分方程
式といいます。
どういう企みか一言でいいますと、初期状態から今までの状態
を元に、現在の結果を少しづつ足して、次の瞬間を予測しよう、と
いったところです。
当然ながら誤差が発生します。これは、未来を予測するのは困
難だ、という現実の世界の教訓とよく合致する意見だと思います。
（１）オイラー法
ここでは工学的に考えやすいように、元の関数を時間 t が経つ
に従って変化する関数 x=f(t,x)として話をします。x'=f'(t,x)の形の
微分方程式(１階常微分方程式)を x0=f(t0)の形の条件(初期条件)
で解きます。
これまでの使ってきた数値積分と同じで、解を求めようとする
区間を微小区間 h に分割し、t1=t0+i・h (i=1,2,...,N)における
x(ti)=xi を順次求めていく。x'=dx/dt を t に関する微係数とする。1
次の微分方程式の一般形は
x'=f(t,x)　初期値：x(t0)=x0 (*1)
です。ti における x(i)が求まれば、式(*1)より x'=f(i)が決まります。
したがって、ti から ti+1 の微小区間 h の積分は
xi+1=xi+h・f(i)　(*2)
で与えられます。
サンプルプログラムは、単振り子の運動です。変位 x、時間 t、
力が復元力 f=-x とした場合、x'=-1・f'(x)と定義される常微分方
程式です。解析解 sin(t)と数値解を比較して表示します。本当は
数値解と解析解の結果を表示させてるのですが、この範囲では
線同士が重なりあって青線しか見えていません。お互いの初期状
態が合致するように注意しときましょう。
この解法には h2・f/2!に相当する誤差が含まれ、時間とともに誤
差（図の黄線）が大きくなっていくことがわかると思います。ちな
みに図では見やすいように、プログラムで誤差を 1.0E+4 倍してあ
ります。誤差を小さくするには h を小さくとる必要があります。それ
は時間刻みを小さくするということで、同じ様な結果を得るにはそ
れだけ計算時間がかかることになります。ちょうどよい時間刻み
を見い出すのがコツとよく言われます。
(GAMBAS)
' Gambas module file
' Euler method (True BASIC => GAMBAS)
' EulerMethodGam.bas
Public Sub Main()
Dim x As Float
Dim v As Float
Dim t As Float
Dim dt As Float
Dim ya As Float
Dim yb As Float
Dim yc As Float
' Define the variables that we'll need
x = 0.0
v = 1.0
dt = 0.01
For t = 0 To 31.4 Step dt
ya = Sin(t)
yb = x
yc = 10 ^ 4 * (Sin(t) - x)
Print t; ", "; ya; ", "; yb; ", "; yc
v = v - x * dt
x = x + v * dt
Next
Stop
End
(Result)
0, 0, 0, 0
0.01, 0.00999983333417, 0.01, -0.00166665833336
0.02, 0.01999866669333, 0.019999, -0.00333306666923
0.03, 0.0299955002025, 0.0299960001, -0.00499897504343
0.04, 0.03998933418663, 0.03999000059999, -0.00666413355839
0.05, 0.04997916927068, 0.04998000209992, -0.00832829241668
0.06, 0.05996400647944, 0.05996500559964, -0.00999120195404
0.07, 0.06994284733753, 0.0699440125988, -0.01165261267283
0.08, 0.07991469396917, 0.0799160251967, -0.01331227527518
0.09, 0.08987854919801, 0.08988004619208, -0.01496994069664
：　：　：　：
31.39 までズラーッと出力されます。
（２）ルンゲ・クッタ法
オイラー法からの続きです。コンセプトは、もうちょっとましな予
測の仕方を考えよう、ということです。
二次のルンゲ・クッタ法は、点（ti,xi)における微係数から x の増
分 ki を推定します。つまり、
k1=h・f(ti,xi)
です。
次に、点(ti,xi)と(ti+1,xi+1)の中間点を考えて、(ti+h/2, xi+ki/2)での
微係数から、これに平行な直線を描いて増分 k2 を推定します。
k2=h・f(ti+h/2, xi+k1/2)
から、二次のルンゲ・クッタ法の式が、
xi+1=xi+k2
で与えられます。k2 の中に、k1 が含まれてるのがポイントです。
四次のルンゲ・クッタ法の場合も同様にして、k1 から k4 までの
量を求めた後、xi+1-xi の誤差が h4 の項まで消えるように重みを
(k1:k2:k3:k4)=(1:2:2:1)として平均した量を k として、x(i+1)=xi+k を
求めます。
xi+1=xi+1/6・(k1+2・k2+2・k3+k4)
で、
k1=h・f(ti, xi)
k2=h・f(ti+h/2, xi+k1/2)
k3=h・f(ti+h/2, xi+k2/2)
k4=h・f(ti+h, xi+k3)
k4 には k3 が、k3 には k2 が、k2 には k1 が含まれてるのがポイント。
解の曲線が連続してるんだったら、それらしい辺りに次の点がく
るはず。という予測のもと、なんとか誤差を小さくしようとしてる努
力が忍ばれます。
今回のサンプルプログラムの結果も、オイラー法による出力結
果と同じ形式で出力します。今回も誤差を 1.0E+4 倍してます。オ
イラー法での誤差の発散具合と比較してください。
(GAMBAS)
' Gambas module file
' RungeKuttagam.bas
' True BASIC => GAMBAS
' 4th Runge-Kutta method (True BASIC)
Public Sub Main()
Dim x As Float
Dim v As Float
Dim t As Float
Dim dt As Float
Dim k1 As Float
Dim m1 As Float
Dim k2 As Float
Dim m2 As Float
Dim k3 As Float
Dim m3 As Float
Dim k4 As Float
Dim m4 As Float
Print "4th Runge-Kutta method m・d2x/dt2=-kx"
x = 0.0
v = 1.0
dt = 0.01
Print "t"; ","; "x=sin(t)"; ","; "x(cal)"; ","; "error"
For t = 0 To 31.4 Step dt
Print t; ","; Sin(t); ","; x; ","; 10 ^ 4 * (Sin(t) - x)
k1 = dt * gFnc(v)
m1 = dt * fFnc(x)
k2 = dt * gFnc(v + m1 / 2)
m2 = dt * fFnc(x + k1 / 2)
k3 = dt * gFnc(v + m2 / 2)
m3 = dt * fFnc(x + k2 / 2)
k4 = dt * gFnc(v + m3)
m4 = dt * fFnc(x + k3)
x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6
v = v + (m1 + 2 * m2 + 2 * m3 + m4) / 6
Next
Stop
End
Function gFnc(v As Float) As Float
Return v
End
Function fFnc(x As Float) As Float
Return - x
End
温度変化を数値積分
ここで元ネタとした論文は、この GAMBAS を紹介するブログ
のテーマ「ロスト・テクノロジーを蘇らせる」によく合致した内容で
す。教育現場で、もとは N88-BASIC で作られた教材プログラムを
VBA に書きなおして再利用しようといった論文です。PDF で配布
されてるので各自検索してダウンロードしてみてください。
しかし、遅いとか読みにくいとか、なにかと評価の低い従来型
の行番号つき BASIC ですが、使ったことのない世代にとやかく
言われる筋合いはないというものです。そこでそれらを「旧制
BASIC」と呼ぶことを提案したい。こう呼べば、バンカラながらも
硬派な印象を受けるのではないだろうか。お勤め先やご家庭で、
ご利用していただきたいと思います。
それにしても「旧制 BASIC」という語感よりも、「旧制
FORTLAN」のほうがより硬派な感じを受けるのは否めませんな。
論文「BASIC リユースと創発」より、
例題 4.6 ルンゲ・クッタ法
File Name: "ex4-8"
均一液相反応の温度変化
引用「よくわかる化学の情報技術」オーム社、83 ページ
このプログラムは常微分方程式の軌道を求めるルンゲ・クッタ
法という数値積分法を使っています。元の微分方程式から少しづ
つ次の点を予測していく仕組みなのでどうしても誤差が積み上げ
られていきます。ここでは、均一液相反応の温度変化を予想する
サンプルを載せます。
反応率 x と温度 T(℃)の間に次のような関係があるとします。
dt/dx=-65.0-15.58・(t-340)/{k(1-x)}
温度 T の初期値は 340℃とします。これは、厳密に解くことので
きない反応率 x の非線型一階微分方程式です。反応係数 k は次
のように与えられます。
k=1.17・1017・exp(-44500/RT)[h-1]
x の範囲が 0 から 0.95 について、x と T の間の関係を計算しま
す。
(GAMBAS)
' Gambas module file
' ex4-8gam.bas
' File Name: "ex4-8"
' 引用「よくわかる化学の情報技術」オーム社、83 ページ
'
' DEF FNK(T)=1.17E+17 * Exp(-22400/(T+273.15))
' DEF FNF(X,T)=-65-15.58*(T-340)/FNK(T)/(T-X)
Public Sub Main()
Dim T As Float
Print "均一液相反応の温度変化"
Print "浴槽温度 T（℃）"
' Print " T = ? ";
' Input T
T = 340.0
RKuta(T)
Stop
End
Sub RKuta(T As Float)
Dim x0 As Float = 0.00
Dim xe As Float = 0.95
Dim dx As Float = 0.05
Dim x As Float
Dim K1 As Float
Dim K2 As Float
Dim K3 As Float
Dim K4 As Float
Print "反応率(-)", "温度変化(℃)"
Print x0, T
For x = x0 To xe Step dx
K1 = dx * FNF(x, T)
K2 = dx * FNF(x + dx * 0.5, T + K1 * 0.5)
K3 = dx * FNF(x + dx * 0.5, T + K2 * 0.5)
K4 = dx * FNF(x + dx, T + K3)
T = T + (K1 + 2 * (K2 + K3) + K4) / 6
Print x + dx, T
Next
End
Function FNK(T As Float) As Float
Dim FNK As Float
FNK = 1.17E+17 * Exp(-22400 / (T + 273.15))
Return FNK
End
Function FNF(x As Float, T As Float) As Float
Dim FNF As Float
FNF = -65 - 15.58 * (T - 340) / FNK(T) / (1 - x)
Return FNF
End
(Result)
均一液相反応の温度変化
浴槽温度 T（℃）
T = ? 340
反応率(-) 温度変化(℃)
0
340
0.05 336.841755390705
0.1 333.920748802926
0.15 331.323483973615
0.2 329.132212160438
0.25 327.405999397189
0.3 326.164620987596
0.35 325.384398882746
0.4 325.008478036539
0.45 324.96508172997
0.5 325.184477066938
0.55 325.609440219045
0.6 326.199250159372
0.65 326.929794658459
0.7 327.792489114165
0.75 328.794080022554
0.8 329.959337582668
0.85 331.340507568081
0.9 333.0465265363
0.95 335.36741332636
参考文献
・吉村忠与志, 青山義弘, 西宮辰明, Excel/VBA を用いた BASIC プログラムのリユ
ースとそれに伴う創発の育成, J.Comput.Chem.Jpn.,Vol.1,No.2, p.59–64, 2002
（３）二階微分方程式
振り子の運動方程式を解き、運動の軌跡を吟味します。一般の
運動方程式は二階微分方程式で与えられます。
それは、
d2x/dt2=f(t,x,dx/dt)
です。ここで、変数 x の時間微分を v と置くことで、二階微分方程
式を二元連立微分方程式の形に置き換えることができます。この
置き換える箇所が今回のポイントです。
dx/dt=g(t,x,v)
dv/dt=f(t,x,v)
という感じです。
この連立一次微分方程式を解けば運動方程式の軌跡を計算
することができます。
四次のルンゲ・クッタ法は次の式で与えられます。
xi+1=xi+1/6・(k1+2・k2+2・k3+k4)
vi+1=vi+1/6・(m1+2・m2+2・m3+m4)
で、
k1=h・g(ti, xi, vi)
m1=h・f(ti, xi, vi)
k2=h・g(ti+h/2, xi+k1/2, vi+m1/2)
m2=h・f(ti+h/2, xi+k1/2, vi+m1/2)
k3=h・g(ti+h/2, xi+k2/2, vi+m2/2)
m3=h・f(ti+h/2, xi+k2/2, vi+m2/2)
k4=h・g(ti+h, xi+k3, vi+m3)
m4=h・f(ti+h, xi+k3, vi+m3)
です。
復元力が非線形-sin(x)の場合と、それを線形近似して-x とした
場合の結果を比較します。
(GAMBAS)
' Gambas module file
' OscillationGam.bas
' True BASIC => GAMBAS
' Linear & Nonlinear Oscillation (True BASIC)
' mode1:-k*sin(x), mode2:-k*x
Public Sub Main()
Dim d0 As Float
Dim x As Float
Dim v As Float
Dim t As Float
Dim dt As Float
Print "Linear & Nonlinear Oscillation f(x)"
Print "initial angle(deg) = ";
Input d0
' 初期値
d0 = d0 * Pi / 180
x = Sin(d0)
v = 0.0
dt = 0.02
' 計算ループへ
Runge(x, v, dt)
Stop
End
Sub Runge(x As Float, v As Float, dt As Float) As Float
Dim t As Float
Dim k1 As Float
Dim m1 As Float
Dim k2 As Float
Dim m2 As Float
Dim k3 As Float
Dim m3 As Float
Dim k4 As Float
Dim m4 As Float
For t = 0 To 28 Step dt
k1 = dt * gFnc(v)
m1 = dt * fFnc(x)
k2 = dt * gFnc(v + m1 / 2)
m2 = dt * fFnc(x + k1 / 2)
k3 = dt * gFnc(v + m2 / 2)
m3 = dt * fFnc(x + k2 / 2)
k4 = dt * gFnc(v + m3)
m4 = dt * fFnc(x + k3)
x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6
v = v + (m1 + 2 * m2 + 2 * m3 + m4) / 6
Print t; ",", x; ",", v
Next
End
Function gFnc(v As Float) As Float
Return v
End
Function fFnc(x As Float) As Float
Dim fx As Float
' mode1
' fx = - Sin(x)
' mode2
fx = - x
Return fx
End
(Result)
時間ステップが同じはずなのに、時間が経つに連れて時間位
相がずれていってしまってるのが見えます。
今回のポイントは変数 x と v の 2 変数の関係を導入したところ
です。これによって単に時系列的な運動の見方の他に、位相面解
析という見方で運動を眺めることができるようになります。そうい
う見方したのが次の図になります。
力学的な運動方程式でいう位相面は普通、横軸に変位 x、縦
軸にその変位 x の時間微分、つまり速度 v の x-v 平面です。
振り子の運動では変位が最大のところでは速度がゼロに、ま
た速度が最大のところでは変位はゼロになります。この軌跡は、
外力もなければ減衰もない、理想状態の振り子がゆらゆら揺れ
ている軌跡です。
時系列プロットと位相面図を比較すると、振り子運動の軌道の
性質は割と一致しているのに、時間が経つ方向に誤差が累積し
ていってるらしいことがわかります。
そこで、次は時間ステップを変えてみたものを比較して見ること
にします。それが次の図です。ここでは初期値 60deg で同じです。
時間ステップが荒いと誤差は大きくなる。けど時間が経つに連
れてというよりも、始めっから誤差がありますな。
一方、初期値をふってみた結果では次のようになります。ここで
縦軸は解析解と数値解の差をとり、時間ステップは dt=0.02 で同
一です。
数値解は小さい θ の場合（sinθ≒θ）のように近似してるため、
振幅が小さいほど誤差の累積が少なくなります。そのとおりの結
果です。
まあ今回の例ではただの振り子を使ったんですが、解析的に
解けない非線形方程式もこの手を使って運動の性質を眺めるこ
とができるという楽しい方法です。
・参考文献
[1] 「シミュレーション物理の功罪」(p.50)
[2] 戸川隼人、UNIX ワークステーションによる科学技術計算ハンドブック、基礎篇 C
言語版、サイエンス社
ルンゲ・クッタ法について先人が工夫したいろいろな方法が挙げられています。
[3] 田辺行人、藤原殻夫、常微分方程式、東京大学基礎工学双書、東京大学出版
会
常微分方程式を体系的に解説してます。FORTRAN77 によるオイラー法、ルンゲ・
クッタ法、ミルン法を計算させるサンプルプログラムが掲載。
あとがき
ブログにアップロードできる１つのファイルの大きさが 1MByte
と制限されているため、PDF ファイルの大きさがそれを超えるまえ
にまとめさせていただきました。
この後に続く非線形方程式の数値解法や解析手法。それとカ
オスやフラクタルについての記事は、後の回にまとめさせていた
だきたいと思います。
「牛と戯れた日々を想って」ブログアドレス
http://acc-holy.cocolog-nifty.com/blog/
このマニュアルの著作権は放棄しませんが、転載は好きなようにしてください。この
マニュアルよりももっとよくできたものが出来るのならそれに越したことはないと思
います。”すでにあるソースは有効に利用する”というのはプログラマの良い伝統で
す。みんなが幸せ。