数値計算法
物理学実験

III(1999)
Fortran 版第１版 (1996.10)
from Ｃ版第１版 (1994.10)
from PASCAL 版改訂第３版 (1990.4.7)& 同改訂第４版 (1992.10.20)
c 富阪幸治 (1990,1992,1994,1996)

目次
この実験を通じて学ぶこと
はじめにプログラム言語について
まず小手調べ（連接）関数と分岐
常微分方程式の数値解法
法
繰り返し計算修正法
ルンゲ−クッタ法
偏微分方程式の数値解法
熱伝導方程式
陽解法
境界条件と初期条件
繰り返しと配列型方程式
連立方程式の解法
ガウスの消去法
数値積分
台形則
無限区間の積分付録
第３節への準備−配列型
グラフィクスエデ ィター（キャンパス情報リテラシ千葉大学情報処理教育研究会編）（サイエンス社）戸川隼人著
注 教科書：ザ・

(0)
(1)
1.1
1.2
1.3
(2)
2.1 Euler
2.2
2.3 Euler
2.4
(3)
3.1
3.2
3.3
3.4
3.5 Poisson
(4)
4.1
(5)
5.1
5.2
( )
1
2
3 muleFortran 77

1

0. この実験を通じて学ぶこと

2

この実験を通じて学ぶこと0

この実験を通じて学ぶことは，３点ある．

 数値計算のいくつかの方法についてそれを身につける．
 プログラムによって，その方法を表現することが出来るようになる．
 データを読む，プログラムを書く，実行する，などの基本的動作をパソコンを使って行え、

Unix やＭＳーＤＯＳなどのオペレーティング・システムについても基本的動作については理解できるようになる．

1

初めにある方程式を解いて答を得ようとする場合に，答が「手」で求められることは，そんなに多くはない．我々のまわりの物理現象は，微分方程式で，表わされることが多いが，この解が，いつも既知の関数で表わされるわけではない．このような時，計算機を用いて解を得ようとする方法
を，数値計算法と呼ぶ．もちろん d2 y=dx2
y の解が yax bx と計算機で求められるのではなく，x なになにの時 y なになにと，数値で求められる，と言う意味である．

=

1.1

=

=0

= sin + cos

プログラム言語について

数値計算法自体は，「 方法」であるから，数学的に記述できるが，実際に計算機に仕事をさせるとなると，計算機にわかる「記述法」で問題を記述してやる必要がある，これを，プログラム言
語という．
ここでは，
という言語を，用いることにする．まず始めに，このについて学びながら， ，等，どの言語でも共通のプログラムの基本的構造を学ぶことにする．

FORTRANC BASIC

1.2

FORTRAN

まず小手調べ（プログラムは順番に）

1.1

例題
a; b ２つの整数をキーボード から入力し，その和，差，積，商を求めて，ディスプレイ
に出力せよ．

ccc

Reidai No.1
integer a,b,wa,sa,seki,shouread(*,*)a,b
write(*,*)a,b
wa=a+bsa=a-bseki=a*b
shou=a/b
write(*,1000)wa,sa,seki,shou
1000 format('wa=',i4,' sa=',i4,' seki=',i4,'

shou=',i4)1. 初めに

3

stop
end
重要点

1. プログラムは，書き出す位置が大事で、一般の文は７カラムめから始まる。始めの５カラム
は文番号を書く欄であり、６カラムめは継続行を示す記号を書く欄である。また最初に c と2.

記された行は、コメント行であり、計算機には何の効果も持たないが、われわれの心覚えのために利用する。
integer は変数の性格を示す．integer は整数，real は実数，charcter は文字型．
変数とは、様々なデータを記憶しておく場所に付いた名前で、電卓のメモリーの名前のようなもの。この例では、の６個の変数を用意し、これには整数を格納することを計算機に教えている。
read，write は読み，書きを行なう．例題では、変数 a と b（入出力並び）に値を読み込み、その値を書き出している。read(*,*) で最初の*は標準に計算機がデータを読む機器であるキーボード から読むことを、write(*,*) の場合はディスプレイに書き出すことを意味する。
=は代入で，右辺を計算して左辺に代入する事を示す．
数学の等号と意味が異なることに注意。
四則は，+，-，*，/を用いる。
（注意）ただし、整数どうしの割り算の結果は切り捨てで整数で与えられる。
変数は、先頭が英字、以降は英数字で 文字以内．
数は、整数は、、 のように、実数はまたはのように書く。
この場合は、二つの整数を「 , 」で区切って並べて入力する。入力するとは、キーボード から打ち込んだあと、リターンキーを押すことである。
様々な名前では大文字と小文字は区別しない。abc と Abc はおなじ名前である。
入出力文で read(*,*) 、write(*,*) の（ ）の２番めの部分はどのような書式（書き方）で、があらかじめ決めた書式を用い
並び a,b を読み書きするかを指定している。*はることを、write(*,1000) はという文番号のついた文で指定された書式で出力することを意味する。
文で'wa=' の部分は wa=という文字を書くことを、そのつぎの、
i4 は並びの対応する変数 wa の値を４桁の整数で書き出すことを意味している。

a,b,wa,sa,seki,shou

3.
4.
5.
6.
7.
8.
9.
10.

123 0

123.4

format1000

1.234E2

FORTRANformat1.1.a 上のプログラムを，ふたつの 実数 の，和，差，積，商を求めるにはどうすればよいか．
( + )2 を計算して、変数 x に代入する文を書け。ただし、a 、b 、c は実数を記憶する問題 1.2
+
問題

b

c

a

a

b

変数であるとする。

重要事項 ここでみたように，プログラムの基本的構造の一つは順番に実行して行く（連接という）という構造である．このほかに後に出て来るが，ある条件によって処理が変わる（分岐 １．
３），おなじ手順を繰り返す（ ２．２）といった基本的構造がある．

x

x

1.3

例題

c

関数と分岐

1.3

exp を求めよ．

実数 x を読み，その２乗，３乗，平方根，サイン，

1. 初めにcc

4

REIDAI 1.3
real x,x2,x3,rx,sx,exread(*,*) xx2=x*xx3=x*x*xrx=sqrt(x)sx=sin(x)ex=exp(x)write(*,*) x2, x3, rx, sx, ex
stop
end

cos(x)，tan(x)，log(x)=log (x)，abs(x) などがある．関数は自分で作ることも注意この他に，できる．

x の３乗を求める関数を作れ．

例題 1.4ccce

REIDAI 1.4
real x,x2,x3,rx,sx,ex
real Cuberead(*,*) xx2=x*xx3=Cube(x)rx=sqrt(x)sx=sin(x)ex=exp(x)write(*,*)x, x2, x3, rx, sx, ex
stop
end

c
c
c

ここからが自分で作った関数

real function
real aCube=a*a*a
returnend

Cube(a)重要事項

1. 関数は x(y) という形をしている。

1. 初めに5

2. x の値が、a に コピー され、その３乗が計算されて、関数名である Cube に代入され持ち帰ら
3.
4.

れる。
もとの主プログラムに戻るには、return を用いる。
プログラム単位（主プログラム、関数、サブプログラム）の終了は end である。

max

例題 1.5a，b の２つの整数のうち大きいものを求める関数
をつくれ．
これには，プログラム構造の分岐を用いる．
で分岐を表現する最も普通の方法は，if 文（教科書）を用いる．

Fortranp.

if (条件) then.... 条件が満足されるときに.... 実行する手順
else.... 条件が満足されないときに.... 実行する手順（これがないときには else 以下は省略）
end if
したがって，この場合はinteger function max(x, y)
integer x,yif (x.gt.y) thenmax = xelse
max = yend ifreturnendとなる．
問題 1.6

二数 a，b を読み込み，その最大値と最小値を出力するプログラムを作成せよ．

問題 1.7

これを三数 a，b，c の最大値と最小値を出力するものに変更せよ．

2. 常微分方程式の数値解法

6

常微分方程式の数値解法

2
2.1

Euler 法

一階の微分方程式 dy=dx「充分」に小さくとると、= f (x; y) を初期条件 x = x0 で y = y0 のもとに解くことにする。h を( + h)  y(x) + hf (x; y);(2:1)
これで、(x ; y ) から (x +1  x + h; y +1 ) の値が求められるので、これを繰り返していけば、
y x

n

n

n

n

n

微分方程式の解が得られる。

y6(x +1; y +1)
#
6
#
#
#
#
#
hf (x ; y )
#
#
##
n

n

n

(x ; y )
n

2.2

n

h

n

-x

繰り返し計算前の節から、微分方程式を解く主要な部分は、

y=y+f(x,y)*hx=x+hであることがわかった。これを繰り返し実行すれば微分方程式の解が得られることになる。x=x+hという文は、「 x に記憶されている中身を h だけ増やして、また x に記憶させる。」という意味であることに注意。
ここで最後の基本的プログラム構造繰り返しがでてきた．
では繰り返しは、Fortrando 文番号 1 i=1,ni の値を 1,2,3,...n と変えて do と
文番号のついた文との間を繰り返し実行する。
文番号 1 continueという文で記述される。i を制御変数と呼ぶ。一つ飛ばしや、二つ飛ばしも

do 文番号 i=1,n,2ようにして書ける。

do 文番号 2 x=x0,xmax,dx

2. 常微分方程式の数値解法

7

x の値を x0,x0+dx,x0+2dx,...xmax と変えて do と
文番号のついた文との間を繰り返し実行する。
文番号 2 continueのように、実数を制御変数に選ぶこともできる。
問題

i=n,n-1,n-2,...,1 のように実行するには、どのように do ループを書けば良いだろうか。

例題 2.1cとして、１から１００までの整数の和を求めてみよう。

Rei 2.1
integer i,wa
wa=0
do 100 i=1,100
wa=wa+i
100 continue
write(*,1000)i-1,wa
1000 format('Sum of 1,2,...',i3,'=',i6)
stop
end

P

=

練習問題
i
i=1;N
ログラムを修正せよ。
練習問題
いか。

P

i

=1

;N

i2

(a)

を使って、この結果が正しいかどうかを確認できるように、プ

(b)

=（

）を計算するには、プログラムをどのように変更すれば良= 0 y （この答えは= 0 = 0 とする。

(1)

課題 2.3dy=dxxけ。初期値は xでy

2.3

Euler 法を用いて解修正 Euler 法

Euler 法の計算法を、
y~ +1 = y + f (x ; y ) 1 hとしたものを修正 Euler 法という。
n

課題 2.4

2.4であるが）、これを

n

n

n

yn+1

課題１をこの方法で解き、（課題= y + f (x ; y ) + f2(x +1; y~ +1) 1 h
n

n

n

n

n

(2:2)

2.3 との）誤差の差を調べよ。

ルンゲ−クッタ法

Euler

この法は、関数の変化が直線的でないと誤差が大きくなる。もう少し誤差の少ない方法
としてルンゲ−クッタ法がある。

(1)yn= y + h2 f (1);
n

n

fn(1)

= f (x ; y );nn

2. 常微分方程式の数値解法
yn(2)
yn(3)
yn+1

2.5

8
= y + h2 f (2);
n

fn(2)

n

= y + hf (3);
n

fn(3)

n

= f (x + h2 ; y(1));
nn= f (x + h2 ; y(2));
nn= y + h6 (f (1) + 2f (2) + 2f (3) + f (4)); f (4) = f (x +1 ; y(3))
n

n

n

n

n

n

n

n

(2:3)

連立微分方程式

連立微分方程式dydx= f (x; y; z)y= f (x; y; z)を解くには、どうしたらよいだろうか？修正 Euler 法を例にとると、
y~ +1 = y + f (x ; y ; z ) 1 h
z~ +1 = z + f (x ; y ; z ) 1 h
f (x ; y ; z ) + f (x +1 ; y~ +1 ; z~ +1 )y +1 = y +1h2f (x ; y ; z ) + f (x +1 ; y~ +1 ; z~ +1 )z +1 = z +
1h
2
dz

n

n

n

n

dx

z

n

n

y

n

n

n

n

n

z

n

n

n

y

n

n

n

y

n

n

n

z

n

n

n

z

n

n

n

(2:4)

(2:5)

とすればよい。プログラムで書けば、do 100 x=x0,xmax,hy1=y+fy(x,y,z)*h
z1=z+fz(x,y,z)*hy2=y+0.5*(fy(x,y,z)+fy(x+h,y1,z1))*hz=z+0.5*(fz(x,y,z)+fz(x+h,y1,z1))*hy=y2
100 continueの様にすればよい。
質問上のプログラムで一度 y2 に新しい y の値を代入した後、２行下であらためて y に代入している。これはなぜか？
連立微分方程式を用いると、２階以上の微分方程式を解くことができる。y 0d2 y=dx2
f x; y; y 0 は次のように書き換えられる。

= (

)

dz
dx

= f (x; y; z)dydx

課題 2.5期条件は x上の置き換えを利用して、y 00
で、y
; y0
とせよ。

=0

=0 =1

= z と置くと

=z

(2:6)

= 0y（答えは言わずと知れた三角関数）を解け。初3. 偏微分方程式

93

偏微分方程式

3.1

準備−配列型

この章では、偏微分方程式たとえば2
 @t@ 9(r; t) = [0 2hm r2 + V ]9(r; t);ihのような式を解くことを考える。この式の微分を差分に変えて解こうとすると、解こうとする場のところの の値が 1 、xのところの の値が 2 、x
所に番号をつけて、xのところの の値が 3 、xのところの の値が 4 、などと表して、その 1 2 3 4のある時刻での値を求めることを考える。この時の準備のためにまずここで，配列型について学び、その使用法について習熟する。
同じ性質のデータが多数ある時，変数名を多数使って ， ， ， ， ， ， ，などとするのは後の取扱に面倒である．こんな時はデータに順番をつけて a1 ; a2 ; a3 ; : : : とする方がよい．このようで取り扱うのが配列型（教科書６６ページ）である．配列型はプログラム上
なデータを
では
とか
というふうに表す（はならをあらわす）．変数が配列型であるこ
とは

9

=1
9 = 41

9

9

9

= 21
9

9

9 = 31
9 ,9 ,9 ,9 ,...,

a b c d e f g

Fortrana(1) a(i)a(i) i=3a(3)integeradimension a(100)realbdimension b(10)a(1) a(2) a(3) ... a(100) 100の様にして宣言する．こうすると，
，
，
， ，
の
個の整数型の記憶場所と，
，
，
， ，
の
個の実数型の記憶場所が準備される．これらやなどの一つ一つを配列要素という。これは、b(1) b(2) b(3) ... b(10) 10
integerreala(1) b(10)a(100)b(10)の様に、宣言しても良い。さらに、負の番号や０をとる配列要素を使うには、

integerreala(0:100)b(-10:10)a a(0),a(1),...,a(100)のように宣言すれば良い。配列 はの１０１個の配列要素が、配列の２１個の配列要素がメモリー上に確保される。

10),b(-9),...,b(10)問題
test.d というファイルに次のようなデータが入っている．

15 クラスの人数
88 １番の人の点数
45 ２番の人の点数
60 ３番の人の点数......
90 １５番目の人の点数

b は b(-

3. 偏微分方程式

10

これを読んで，平均値，標準偏差を求め，元のデータと共に出力するプログラムを作れ．
注意１：平均値，標準偏差の求め方は教科書６６ページ参照
注意２：キーボード からでなくというファイルからデータを読むには，（ 教科書１７９ページ）つぎのようにする．

test.dopen(2,file='test.d',status='old')
read 文の前にこれを実行する。
意味：test.d という既に存在するファイルからデータを読んだり、そこにデータを書いたりすることを示す。
status='new' とすると、新しいファイルをつくってそこへ書き込むことを示す。

...
read(2,*) a(i)実際に読み出し，a(i) に記憶する．
read(*, を read(2, とすると、キーボード から読み込む代りに open 文で対応づけられたファイル test.d から読む。必要なら format の指定を行なう。
書き出す場合は write(2,... のようになる．
問題 前の問題ができたら，さらに，クラスの最高点と最低点を求めて元のデータと共に出力するプログラムを作れ．
参考：最大値の捜し方

integer saidai,i,n...
saidai=0
do 100 i=1,nsaidai=max(saidai,a(i))
またはif (a(i).???.saidai) thensaidai=???
else??????????
end if
100 continue
質問

この？には何が入るか。

3. 偏微分方程式
3.211

熱伝導方程式

この節では、熱が高温の物体から低温の物体へながれる有り様を表す、熱伝導方程式を取り上げ、偏微分方程式の差分法による解き方を学習する。鉄の棒の片方の端を，突然高温に熱したとする。熱伝導で，その端から他方へ，熱が流れ，しだいに全体が高温になるだろう．このときの温度の変化は，場所場所でどの様に変化するだろうか．これを記述する方程式が熱伝導方程式である．

図のように，棒の長さの方向に，x 軸をとる．単位面積の断面積を通って流れる，熱エネルギーの量 q は，熱伝導率  と温度勾配 @T =@x の積で表される．

()

(

)

q

(3:1)

= 0 @T@x

15  50Js01 m01 K01 である．単位体積当りの熱エネルギー E の時間

ちなみに鉄の熱伝導率は，
変化は

@q
=
0
@t
@x

(3:2)

@E

で与えられる。質問：何故か？
ここで、棒の物質の単位体積あたりの比熱を c とすると，E
温度に依らずに一定とすると，温度変化の式は，
c( ) =  @2 T (x; t)@t@x2@T x; t()= cT だから，比熱と熱伝導率が
(3:3)

()

となる．これが、熱伝導によってある場所 x の温度 T x がどのように時間変化してゆくかを表す熱伝導方程式である
さらに、

t0

(3.3) 式は，

と置くと，

x0

= Lx

u0

= TT

0

= (cLt2=)

@u0

2 0

@ u
=
0
@t
@x02

(3:4)

3. 偏微分方程式12となる．x0 、u0 、t0 はすべて、次元のない量になっている。このように、方程式を次元のない量（無次元量）に対するものに書き換えることを無次元化あるいは規格化という。
は企画化された熱伝導方程式である。

(3.4)

質問 棒が上に凸の温度分布をしていた時、棒の外から（の外へ）熱の流入が（流出）ないとす
ると、中心の温度は下がるか上がるか。これは方程式のどこから理解できるか。また下に凸の温度分布をしていた場合はどうなるか。

3.3

陽解法

(3.4)
()
+1

書くのに面倒なので，これ以降
式の 0 は省略する．これを，この節では，差分法で解く．
x の区間x上に等間隔 h に N
個の点を考え，そこでの温度変化を追跡する．温度を u0 ; u1 ; u2 ; : : : ; uN と書く．点の間の間隔が小さければ，微分を差分に書き直せる．つまり，

0 1

@u
@x
@2u

@x2

= u +1h0 u
i

(3:5)

i

= u +1 0 2hu2 + u 01
i

i

(3:6)

i

よって

(3:7)
( + 1t) = u (t) + 1h2t [u +1 (t) 0 2u (t) + u 01(t)]右辺は時間 t の時の値のみで計算できるので，1t 後の x での温度が計算できる．これを，0; 1; 2; : : : ; Nについてやれば，1t 後の温度分布がわかる．
ui t

i

i

i

i

i

3.3.1

境界条件と初期条件

0 C= 0 for = 1 2
100 C( = 100)i; ; : : : ; N である．これを，初期条件 と呼最初，この鉄の棒は  であった．uiぶ．また，片方の端を加熱した後，高温の端は，
低温の端は，  に保たれていたとすで uNである T0．これを 境界条件 という．
る．いつも u0

=0

3.3.2

=1

0 C

繰り返しと配列型

u0 ; u1 ; : : : ; uN のような，変数の集まりをFortran で使うには，配列 を用いるということは既に学んだ．配列は，u(1) とか u(i) といったように使うが，最初に，個数（何番まで使うか）と型（実数か整数か）を指定する必要がある．

real u(0:10)
これは，０から１０まで、１１個の実数の配列要素を使うことを示す（復習）．
のような繰り返しは，
さて、式

(3.7)

do 200 i=1,nu(i)=u(i)+r*(u(i+1)-2*u(i)+u(i-1))
200 continue

3. 偏微分方程式

13
=0

=

の様に書くことができる。iと i N では，いつも ui は決まっているから計算する必要は無2い．ここで r は t=h である．
質問 この上に示した，プログラムには，計算しても正しくできない（間違った答えが出る）点がある．それはなにか？i を１，２，３，と変えながら，どういう計算が行われるか検討してみよ．
課題
この問題を，プログラムを作って解け．t: ， : ， : の時の，温度分布を示せ．rの値はどの程度にとればよいか？

1

3.2

3.4

=01 02 05

ポアッソン方程式

重力や，静電気力などの逆２条に比例する力を及ぼすポテンシャルは

2

2

2

div grad  @@x2 + @@y2 + @@z2 = 4G

(3:8)

に従う（静電気力の時は右辺は =0 になる）．
また，熱伝導の問題で，熱源からの発熱 H を含めると，
c

2

@ T
=
 2 +H
@t
@x

@T

(3:9)

という方程式が得られるが，この発熱を含む熱伝導の定常解（時間的に温度分布の変わらない状）は，
態：@u=@t

0

@2u
@x2

= 0H

(3:10)（ここで，H はもとの H を T0 =L2 で割ったもの）という方程式を解けばよいことになる．これ
式の独立変数が x のみの場合にあたっている．
は，微分を差分に書き直すと，この式は

(3.8)

ui+1

=1

0 2u + u 01 = 0H h2
i

=0

i

(3:11)iとなる．ここで境界条件，u0，uNを付加すると，Hi は与えられているので，この式は連立方程式
u1 ; u2 ; : : : ; uN 01 に対して N
本の式があるとして，解くことができる．

||

002 1 0B
1 02 1B
B
B
0. 1. 02
B
B.. .. . . .
B
B
. .
B
@ .. ..

001||

0 1 0 u1 1 0 H1 h2 + u0 1
0 : : : 0 CCC BBB u2 CCC BBB H2 h2 CCC
1 : : : 0 C B u C B H3 h2 CCC...
. . . . . . ... CCC BBB ...3 CCC = 0 BBB
C
C
B
C
B
C
...
. . . . . . ... CA B@ ... CA B@
CA
2
0 : : : 0 1 02 u 01
H 01 h + u
:::

:::

N

N

(3:12)

N

このような方程式には、初期条件は課されず、境界条件のみが課される。このような方程式を
楕円型方程式 と呼んでいる。この方程式（ポアッソン方程式）を解くために次の節で，連立方程式の解法について学ぶことにする．

4. 連立方程式の解法

14

連立方程式の解法

4
4.1

ガウスの消去法

解くべき連立方程式が

=

Ax

0
a11
B
B a21
A=B
.
B
@ ..

であるとする．ここで，

an1

(4:1)

b

a12

:::

a22

:::

an2

:::

... . . .

0 1
x1
BB x2 C
C
C
x = B
.
B
.
@ . C
A

1
C
a2 C
... CCA
a1n

(4:2)

n

ann

(4:3)

xn

0 1
b1
B
C
B
b2 C
B
b = B . C
@ .. C
A

(4:4)

bn

8 (1)
(1)
(1)
(1) #1>
a11 x1 + a12 x2 + 1 1 1 + a1 x = b1
>
>
(1)
(1)
(1)
(1)
>
< a21 x1 + a22 x2 + 1 1 1 + a2 x = b2 #2
(1)
(1)
(1)
(1) #3a31 x1 + a32 x2 + 1 1 1 + a3 x = b3>>>
: : :(1):(1)(1)(1) #na 1 x1 + a 2 x2 + 1 1 1 + a x = b
消去の第１段では，（ ４．５）式の２番目以降の式から x1 を消す．まず#1 に
a21
m21 =
a11
である．

n

n

n

n

n

n

n

n

nn

n

(4:5)

n

#2 から引くと，
(1)
(1)
(1)
(1)
(1)
(a(1)
22 0 m21a12 )x2 + 1 11 + (a2 0 m21 a1 )x = b2 0 m21b1となる．同じように，#1 に
a31
m31 =
a11
を掛けて#3 から引くと，
(1)
(1)
(1)
(1)
(1)
(a(1)
32 0 m31a12 )x2 + 1 11 + (a3 0 m31 a1 )x = b3 0 m31b1

(4:6)

を掛けて

（４．５）式は

n

n

n

n

n

n

8 (1)
(1)
(1)
(1)
>
a11 x1 + a12 x2 + 1 1 1 + a1 x = b1
>
(2)
(2)
(2)
>
a22 x2 + 1 1 1 + a2 x = b2
>
<
(2)
(2)
(2)
a32 x2 + 1 1 1 + a3 x = b3
>
>
:::
>
>
:
(2)
(2)
(2)
a 2 x2 + 1 1 1 + a x = b
n

n

n

n

n

n

n

nn

n

n

#10
#20
#30 2
#n0

(4:7)

4. 連立方程式の解法(2)15(2) は，次式で定義される．

ここで，aij と bi

( (2) (1)(1)
a = a 0 m 1 a1
(2) (1)(1)b = b 0 m 1 b1

#3'

ij

ij

i

i

i

i

= 2; 3; : : : ; ni = 2; 3; : : : ; n

(4:9)

i; j

j

次の段階では，以降の式から x2 を含む項を消去する．以下同様にして，次々に xi を消していく．
この消去の手順をまとめると，

do 300 k=1, n-1
d=a(k,k)do 200 i=k+1, nm(i,k)=a(i,k)/ddo 100 j=k+1, na(i,j)=a(i,j)-m(i,k)*a(k,j)
100continueb(i)=b(i)-m(i,k)*b(k)200 continue
300 continue
この操作で，元の方程式は0 a(1)
11
B
B
U =B
B
@

Ux

=

(1)
a12
(2)
a22

b

:::
:::

...

(4:10)

(1) 1
a1
(2) C
a2 C
... CCA()an(4:10)nnnn

0 b(1) 1
1(2)
B
C
b
B
C
2
B
b = B . C
.
@ . C
A
(
)
b

(4:11)nnという，形に変換されたことになる．U は上三角行列である．ここまでの操作を，前進消去 という．
（４．１０）式を解くことはたやすい．つまり，n 行目の式から，
xn

これを，n

= b( )=a( )nn(4:13)

nnn

0 1 行目の式に代入すると，
xn01

この二つを n= (b( 0011) 0 a( 0011) x )=a( 0011) 01
n

n

n

n

(4:14)

n

n

n

n

n

0 2 行目の式に代入すると，
xn02= (b( 0022) 0 a( 0022) 01x 01 0 a( 0022) x )=a( 0022) 02
n

n

n

n

n

n

n

n

n

n

n

n

以下同様で，下から順番に xi を求めることができる．まとめると，

n(4:15)

4. 連立方程式の解法

16x(n)=b(n)/a(n,n)
do 200 i=n-1, 1, -1x(i)=b(i)do 100 j=i+1, nx(i)=x(i)-a(i,j)*x(j)100 continuex(i)=x(i)/a(i,i)
200 continue
この部分を，後退代入 という．
課題５ 連立方程式を解くプログラムをつくり，それを利用して（３．１２）式を解け．但し，熱源の分布としては
x
では H
とせよ．

0 1

=1

5. 数値積分

17

5

数値積分

5.1

台形則

定積分
I

=

Z

b

()

(5:1)

f x dx;aX

は= 1lim!0 1xf (x )(5:2)
で与えられる. x = a 、x = b 、x 軸と y = f (x) で囲まれた面積を求めれば良いわけだが、もっとも、簡単には [a; b] を N 等分して，(h = (b 0 a)=N )I ' 1x (f (a) + f (a + 1x) + f (a + 21x)+; :::; +f (b 0 1x))
(5:3)
I

あるいは
I

i

x

i

' 1x (f (a + 1x) + f (a + 21x)+; :::; +f (b))(5:30 )のように求めれば良い。
質問：

これはどのような面積を求めていることになるか。グラフに示せ。

ここでは、もう少し精度の高い「台形則」と呼ばれる方法を考える。定積分I[ ]を求めるのに， a; b を N 等分して，
h=Z

f x dx;()(5:4)= (b 0 a)=N

(5:5)

a

I を次のように近似するのが，台形則である．

"

b

#

1 f (a) + X01 f (a + nh) + 1 f (b) :
I =h
2
2
=1
N

N

(5:6)

n

質問： 「台形則」という名前を参考にして、これはどのような面積を求めていることになるかを考えよ。また、それをグラフに示せ。

はじめから，N を決めて計算するのももちろんよいが，しだいに，N を増やして行って適当な精度で積分値が求められたら，終わるというプログラムが望ましい．

5. 数値積分18

5.1 = 1

= 2, = 4,

課題
N
から N
N
: : : と分点の数を２倍に増やして行き，IN と I2N の差が （これは十分に小さい数）よりも小さくなったら，I がもとまったとして結果を出力するプログラムを作れ．これを用いて、次の積分値を求めよ．
I

Z1= 01 1 +2 x2 dx(= )(5:7)

進んだ課題：余裕のある人向け上の問題では，前に計算された IN の値を使わないで次の I2N を計算した．これは，効果的ではないので，前の値を使うように改良することを考える．
これには，次のようにやればよい．N
のときは，

=1

= [f (a) + f (b)] b 02 a

(5:8)


X 
1
= h f a + ( n 0 2 )h
=1

(5:9)

I1
N を倍にして(h を半分にして)I2 を計算するには，
N

N

JN

n

h

= (b 0 a)=Nを求めて，IN と平均すればよい．つまり，
I 2N= I +2 JN(5:10)

Nj 0jとすれば，刻み幅が，半分の台形則 I2N になる． IN I2N <  という条件を満足したら終わりと
06 とか 小さい正の数である．これを用いて、課題いうことにすれば，良い． は適当なをより効果的なプログラムに書き直せ。

(10

5.2

)

5.1

無限区間の積分
I

()

0

Z1

= 0 f (x)dx(5:11)

で，f x が非常に早く に近付く様なものを考える．

( )=h

I h

1
X
f (nh)=0(5:12)n( )こんどは，下から順番に加えて行って関数の値 f nh が，それまでの，部分和の  倍にまで減っ

sin exp(( + 1) ), (( + 1) )()x

2 のように， の大きい所で，零点を持つものもあるのたところで打ち切ればよい1 。 x
h f n
h も充分に小さくなっていることを確かめた方がで，f nh だけでなく，f n
よい．この手続きで，I h ある h に対する値が求められる．
前の問題と同じように，h を，半分にしていって， I h
I h=が充分に小さくなったら解が求まったと考えよう．
課題 5.2 次の積分値を求めよ．

( )

1

0 x2

j ( ) 0 ( 2)j

ここでは幅１高さ f (nh) の面積と部分和で表される面積を比べている．

5. 数値積分

19
I

Z1

2

p

= 2 0 exp 02x dx(= 2)

この章の参考文献：数値計算の常識，伊理正夫，藤野和建，共立出版．

(5:13)

5. 数値積分20

A. グラフィックスの使い方

Unix

telnetパソコンから、マシンにを用いてリモートログインしている時には、画面には文字しか出力することができない。
グラフィックス（図など）は 端末（画面）に出力する。物理測定機室でワークステーション（、）にログインすると ウインドウシステムが画面にあらわれる。したがって、もっとも簡単にグラフィックスを画面に出す方法は、物理測定機室でワークステーションにログ イン
することである。実習室のウインド ウズパソコンの１台には、パソコンを 端末（ 画面）にするためのソフトがインストールされており、そこに出力することも可能である。

quasar galaxy

XX(edpc20)X

A.1 グラフィックスのプログラム

SuperMongo

まず，図を書く方法にはいくつかあるが、ここでは（スーパモンゴ）を使う。これは、のプログラムのなかで、が提供する関数副プログラムやサブルーチン
副プログラムを呼ぶことによって実現される。例えば座標軸を書くための関数副プログラムなどがすでにシステムで用意されている。
まず、始めに、

FortranSuperMongo

iopen=sm_device('X11 -title tomisaka -g 750x750+100+100')
call sm_graphics()
call sm_defvar('TeX_strings','1')
call sm_erase() sm_device はどこに図を書くかを指定している。この場合は、X ウインド ウ画面上に出力する。

 -title tomisaka の部分は画面上のグラフにタイトルをつけているもので、随時変更すれば良い。

 このとき、X ウインド ウ画面は今自分が座っているワークステーションの画面に出力される。これを、例えば、実習室のウインド ウズ・パソコンに出力する場合は、以下のコマンドをした後で実行すること。

login

setenv DISPLAY edpc20.ed.niigata-u.ac.jp:0.0
これによって、実習室の文字しか出力することができない端末からも、実習室のウインド ウ
ズ・パソコンに出力することによって、グラフィックスを見ることができる。
画面でなく、プ リンターへの出力などは次のように指定する。

iopen=sm_device('postscript')
iopen=sm_device('postfile test.eps')(

上の場合は測定器室のプリンターへの出力、下の場合は図をポストスクリプト形式のファイル 名前がにして保存しておくことを指定している。ただし同時に２つ以上の出力機器を使うことはできない。これ以外の３つのサブルーチン副プログラムの呼び出しは、それぞれsm_graphics：

test.eps)

5. 数値積分21

グラフィックスを書くモード への変更、sm_erase：画面の消去などであるが、そのまま書いておいて欲しい。
さて、これからグラフを書く段階にはいる。

call sm_limits(0.0,1.0, -0.1,1.1)
call sm_ticksize(0.1,0.2, 0.05,0.1)
call sm_xlabel('\\it x')
call sm_ylabel('\\it u')call sm_box(1,2,0,0)
これらの手続きを説明する。

 call sm_limits(0.0,1.0, -0.1,1.1) でグラフを書く範囲が 0 < x < 1 と 00:1 < y < 1:1の範囲に限定される。それ以外の部分は点があっても描かれない。

 call sm_ticksize(0.1,0.2, 0.05,0.1) で、x 軸には 0.1 きざみで目盛をうち、0.2 きざy

みで数値を書き、 軸にはている。

0.05 きざみで目盛をうち、0.1 きざみで数値を書くことを指定し call sm_xlabel('\\it x') は x 軸に x と物理量を書き、
 call sm_ylabel('\\it u') で y 軸に u とその物理量を書いている。
 最後にcall sm_box(1,2,0,0) で座標軸を実際に書いている。引数の最初の１は x 軸の物
理量 x を縦向きに書くこと、次の引数の２はている。ここまでで、枠が書ける。

y 軸の物理量 u を横向きに書くことを指定しA.1.1 線グラフを書く

線グラフはsm_conn を用いて書く。

call sm_conn(x,u,N)xyN(x(1),u(1))
N

ここで x，u はそれぞれ、グラフの点の 座標、 座標の値が入っている配列で、、、 、、がプロットされる 個の点の座標に当たる。最後の は点の数である。
（ただし、 や が０番めから使われている時など は最後の点がであれば、点の数は ではなくであることに注意せよ。）(x(2),u(2)) ... (x(N),u(N))
x ux(0),u(0)x(N),u(N)N
N+1[real x(0:N)]

A.1.2 終了処理

call sm_alphacall sm_redraw(0)X ウインド ウに図を描いている時には、call sm_alpha でグラフィックスを書いていたモード から抜けるとともに、call sm_redraw(0) を呼ぶと終了処理ができる。端末でリターンキーを叩くと画面が消える。

5. 数値積分22

call sm_hardcopy
ポストスクリプトプリンターに出力している時は、call sm_redraw(0) の代りに、call sm_hardcopy
を呼ぶと、プ リンターにグラフが出力される。

A.2 翻訳と実行

スーパーモンゴを使ったグラフを書こうとするプログラムは

f77 プログラムの入っているファイルの名前 -lplotsub -ldevices -lutils -lX11のように、スーパーモンゴのサブルーチン副プログラムが入っているライブラリから必要なサブ
ルーチン副プログラムをとってきて組み込む（ リンクという）作業をすることを指定する。 の後ろの名前がそれぞれライブラリ名である。
これは、~student/.cshrc にすでに設定してあるので、-lcd
cp ~student/.cshrc .cshrc
source /.cshrcというコマンド を一回実行しておけば、以降は、「だけで翻訳することができる。

f77 プログラムの入っているファイルの名前」

B. 環境の整備とレポート の提出方法

Unix

一回だけ
上の環境を作るために以下の作業を行なって下さい（以前からを使っていた人はしなくて良い）。次回からは必要ありません。

quasar や galaxy

cd ~itoh/HOME0./setup.sh -mh
続いて以下を実行して下さい。

cd
cp ~student/.cshrc .cshrc
source /.cshrc
これも次回からは必要ありません。
コマンド で行ないます。プログラムがさて、レポートの提出は、ルに格納されていると、

teishutu

mon1.f というファイ

teishutu mon1.fとすると、プログラムの実行が始まります。そしてプログラムと計算の結果は、画面に表示された後、電子メールで所定のアド レスに配送されます。このとき、データを入力する必要がある場
合は、画面に数値などを入力しエンターキーを押して下さい。次回の実習までには返送されてくる予定です。