数値解析における授業方法改善の試み

(1)

NDC 418， 375．12

数値解析における授業方法改善の試み

宮地 ^潜＊

（昭和54年5月2日受理）

Improvement of Teaching Method in Numerical Analysis

Isao MlyAJI

（Received May 2， 1979）

In the instruction of numerical analysis， this paper proposes a teaching strategy by which students become skillful

in 吹f rogramming necessary for using a computer and in choesing numerical methods suitable for each problem． The instruction in numerical analysis is given as follows． First each numerical methed is taught． Next the FORTRAN program for the method is taught． Finally the practi．ce in preparing the program is performed． The coefficients of functions or input data in the program are changed by each student． Such an educational technique is seemed to be more effective than the instruction by lecturing only．

1．序

論

情報化社会に対応できる幅広い層を養成するために，全ての教科に情報処理との関連に応じて清報処理教育を取り入れる必要があると言われて久しい。電子計算機が設置されさえずれば，情報処理教育ができるような幻想が抱かれている。しかし，どのような道具でも存在するだけでは何の役にも立たない。その道具をどのように使うかは，それを生きたものにするかどうかを左右する。電子計算機についても同様であって，電子計算機をどのように利用するかは，重要な問題である。その意味で，暗中模索の情報処理教育については，多くの試行の繰り返しがより良い方法を見出す近道かも知れない。著者はそのような報告をいくつかした1）〜4）。更に，情報処理教育を前進させたいという願いからこの報告は企画されたものである。一方，高専では学習者の学力の多様化などにより，今まで通りの授業方法では，対応できなくなっている。この面に電子計算機の利用をそれぞれの教科に応じて工夫することによって，質の良い情報を得て学習効果をあげる道があるように思う5）。

ここでは，著者の担当する数値解析の授業を工夫し，一応の成果をあげていることを報告し，情報処理教育および教育方法改善の参考に供したい。

＊機械工学科

プログラミングの好適題材として数値解析がある。しかし，本校では従来数値解析はプログラミングに余り利用さ

れていない。プログラミングはFORTRAN入門だけで終っ

てしまって，数値解析はもちろんのこと，他の教科でもほとんどふれられていなかったようである。プログラミングの知識と数値計算：法の知識を有機的に組合せることによって，これから出合ういろいろな問題に最も適切な計算方法を行なうことができるようになる。そこで，数値解析の授

業において，電子計算機を利用する上で必要なプmグラミ

ングおよびそれぞれの問題に合った計算方法に習熟するように考えた。すなわち，数値解析の各種の手法を説明する

だけでなく，その手法をFORTRANプログラム化して説明

し，プログラムを作成させる演習を行なうように授業を展

開している。例題で用いている関数データなどは各学習

者の固有の学籍番号，および出席番号などを用いて変えて

演習させている6この方法は講義だけによる受身的な授業

に比し積極的に演習する姿が多くみられ，効果がかなりあるように思われる。

以下では，まず授業方法の概略を説明する。次に，数値

解析の授業で取り上げている内容，．すなわち，方程式の解

法，連立一次方程式の解法，最小二乗法，補間法，数値積

分法，常微分方程式，およびモンテカルロ法の代表的解法

を1．つずつ取り上げどのように演習を行なっているかを順

(2)

Table 1 Summary of the conte．nts of the methods in numerical analysis．

Solution

of

equation

Simultaneous linear equation

Least square

method

Interpolation

Numerical integration

Ordinary

differential

equation

Monte Carlo

Method

． Simulation

Beirstow s method

＊ Newton−Raphson s method

von Mises method

＊ Regula falsi method

Bisection method

M・・…P・・at・・n｛1鵠尊詠Mpy）

inverse matrix

＊ Gauss−Jordan s method

Gauss−Seidel s method

翻黙黙翫）

Lagrange s method Newton s method Aitken s method Aitken−Neville s method

．＊ Trapezoidal formula

＊ Simpson s 1／3 formula

＊ Simpson s 3／8 formula

Gauss method

＊ Euler s method

Modified Euler s method Euler−Gauss method

Rtinge−Kutta smethod（（（（ilggco：i， dtC／i

Runge−Kutta−Gill s method

＊Milne s method

＊ Generation of random numbers

・n・・g・al｛（Circle）（Multiple），

Buffon s needle

＊ Simulation for mechanical vibration

No． of pro−

grams

23322君0噌⊥−占−占−−−占211

11⊥11 444一4

444444453242

2

No． of

Steps

29022834449440187576 9210586663343644

47．

49 52 134 52 T9 T5 T4 T5 T6 U9 ﾐ

_63215354

117 Data linage

3229臼ウ臼

7744433442655888966666666111占−

4 No． of presenters

75一75

8912044575 7・964凸 3127．77

69 6FO4 K＝一46扉

ん一Qり1

62 Rate of completion

（％）

87︸84

90 W4 X6 W9 W2

899998

器＝﹁93︸96

1289932

FひFOO

Average

No． of

jobs

3．4

2．7

30902ウ一9臼19一2 Qりρリ

ロ

コ

32

420

ロコ222

弓＝至一β 22． 1

2

212 081

2．3

次説明する。最後に，演習結果のまとめを示し考察を加え．

る。

2．授業方法の概説

この授業は，機械工学科4年生の2クラスについて，・

「数値解析」という100分の後期の授業であるび数値解法の説明をする内容の概要をTable 1に示す。 Fig．．1の演習手順に従って授業方法を説明する。授業の前半に各数値解法の説明を行なうと共にそれぞれのプログラム，入力デ

ータおよび処理結果を説明する。特にプinグラAで工夫す

べき点などは強調する。後半は，説明した解法の1つについてプnグラムを作成させる。プログラムは、関数副プログラムにしてある関数の係数．入力データなどを出席番号または学籍番号で置き換えて作成する。従って，．．得られる数値は，学習者全員異なる。また，プログラムは原則とし

てFig．2のようなFORTRAN形式のマークカードを用いて

カード．ディックを編成することにしている。ただし，Fig．

3のようなパンチカードで編成を希望する者には，．演習の日までにカードをパンチしておくことを条件に許可している。コントm一ルカードなどを付けて，カードディックを

巳@ Start

@P「esentation of

≠氏@alterable Explanatioロ of a metbod

@ and the rogram

@ part

Construction of

≠垂窒盾№窒≠香@deck

Request of 狽??@computation

Modification

盾?@mark car4S

「 Produ．ction of

≠モ??モ求@card サ

kP碧2鯉畦艘j t

m】丞菱u亟：ユ．

@ Is the

@ result

@ right ワ

駈odification 盾?@the program

no

yes Producti6n

盾?@a report

〈調。譲9・ ^no

？

yes Grading Stop

Fig．1 Procedure for the instruction in numerical

analysis．

(3)

数値解析における授業方法改善の試み宮地

Fig．2 Mark card for preparing FORTRAN programs．

、響 F。RT剛S．＾了。M控目＿

閑 1 、

iiiii灘：；1灘ll瓢：蕊：1瓢：1；1灘灘鱗1蹴ll

iiiil：iiiiiiiiiiiiiliiiliiiiiiiiiliiiliiiliiliiiiiiiiiliiiliill灘i

：1・llllFlllltlll／iillllilllllllllllilllll／III：，；ll；・ lll；iSlll＄・ lll￥・ III：i III：i ：1；，llllE

l：：：：，：：：：1：：：1：：：：定：：：：1：：：1：：：：鶉：：：翼：：：霊：：：罵：：：宝：：1罵：：：竃：：：霊：：1：：：：腫

● 鵬 Fig．3 Punch card for preparing FORTRAN programs．

編成して，電子計算機室の受付にて計算依頼する。オペレータによって処理された結果は，返却棚にカードディックと共に返却される。学習者は，その返却された結果を受取り，正しい結果が得られておれば，Fig．4のようなチェッ

Fig．4 Mark card for checking practical results．

クカードの該当する個所をマークし，得られた結果の印刷されたしP用紙の余白にFig．5のようなゴム印を押して必要事項を記入し，カードディックと共に報告書として提出する。教師は，それらの報告が良ければ受取るが，悪ければ再提出を要求する。もし結果に誤りを含んでいる場合，

プログラムを修正すると同時にLP用紙とチェックカードを提出し，再び計算依頼することになる。

課鵬号感激番雛回

年組番氏名

考察

理饒されたことﾃいたこと

Jグ弘の比較エラーの膿因

ﾈど

Fig．5 Rubber stamp required for entering in the program list．

また，数値解法の区切りごとに，精度，収束回数，などについて各解法を比較し，解説している。そして，数値計算における注意すべきことを説明するように努めている。

3．演習方法

数値解析の授業で扱っている内容は，次の8つに大きく分けることができる。（1）方程式の解法，｛2漣立一次方程式

の解法，（3）最小二乗法，〔4）補間法，｛5）数値積分法，（6）常微

分方程式の数値解法，（7）モンテカルロ法，⑧シミュレーシ

ョン。この節では，これらの内容の中から，それぞれ代表的な解法を1つずつ取り上げて演習方法を少し詳しく説明

する。

3．1方程式の解法

技術関係の問題を解くとき，方程式の根を求めることがよくある。方程式の解法には，Cardano法， Ferrari法，

Brown法， Newton−Rapthscm法， von Mises法，はさみうち法，二分法，Horner法， Bairstow法， Lin法， Lehmer 法，Bernoulli法， Graeffe法，逐次近似法など多くの方法がある。中でもNewton−Raphson法は原理も簡単で，プログラムも短いステップでよい。その上，この解法は高次関数，超越関数および複素関数のような広範囲な方程式の根を求めるのに適用できるすぐれた特徴を持っている。このような点から，方程式の解法のプログラミングには，

Newt（，n−Raphs㎝法を取り上げる。

New㎞一Raphscm法によるノ（X）＝Oの根を求める漸化式は，Taylor級数のxに関して1次までの項によって求めら

れ，

xi．i＝xi一一ff；g（ tfiii）lr，（i＝o， 1， 2， …｝・・）

で表わされる。ここで，ノ ④は関数f（X）の導関数である。この式を用いて，適当な初期値Xoを与えることによって，x1， x2，……，Xnを求める。この反復計算の停止条件としては，1つ前に得られた値と新しく得られた値との相対誤差がある値εより小さくなることを用いる。すなわち，

1（Xn＋1−Xn）／Xn＋11＜εである。これらに基づいて， Fig．

6のようなプログラムが作成できる。第2〜10行は，表題

の印刷および初期値の読み込みなどの準備部分である。第 11〜21行が反復計算して根を求める部分で，1回計算することに得られた近似値を印刷するようにしてある。第22〜

27行は，求められた根を印刷する部分である。第30〜34行は，根を求める関数を定義した部分であり，第35〜39行は，その導関数を定義した部分である。

根は，超越関数f（x）＝・ex2−beXについて求める。この関数において，a＝（学籍番号の下2桁）およびδ＝（出席番号）

として各学習者ごとに異なった関数の根を求める。初期値

は，例題と同様にXO ＝・ 1とする。また，停止条件は，ε＝

(4)

ユ234567890123456789U五23ら567890工23456789

ユー−三ーユー11122Z2Z222223333333333 C 砦曽誓 NEWTON−RAPHSON METHOD WR！Tε曜49600G亘

6000 FORMATtlH ，，一 NEWTON一，RAPHSON HOVt，

蓑 1 軸・〃｝

C ．．． 5HOKr．CHr NO YOMTKeMl iO READt3esooO） xa

5000 FORMAT｛FIO．O）

IFtXO．EQ．999q．q｝ STOP

WRITEC4，6e20）

6020 FoRMATclH e， Nn，tgxt−xt−t2XttFtXle）

C ．．． HANPUKU KErsAN

DO 120 1＝1，50

FX累F XO，

NO＝：口1

WRITEC4，6100｝ NO，XO，FX 6100 FORMATCIH ，15，IP2E15．6｝

X冨XO−FX！DF XO】

IF｛AB5こ X一×0，！X⊃・しE●1●OE卿6，

砦 GO TO 200

xo＝x 12n CONTrNUE

C ．．． KON NO INSATSU 200 FX＝FtX）

wRrTEt4，6200） r，x，Fx

6200 FORMATCIIIH ，，HANPUKU KArsU ＝，，ISI

一 IH t−X ＝一elPEIS．6t

畳 1H ●撃F【X，冒・，E工506！！｝

GO TO 10

END

C砦穏蒼KANsU NO 7EIG：

FUNCTION F｛X｝

F竃2000曇X管黄2r1000曽EXP｛X，．

RETURN

ENO．

C ee−f DOUKANSU NO TEIGT

FVNCTION DF（Xl

DF＝4e．o−x−le．o−ExpCx）

RETURN END

一 NEwTON．RAPHSON OU k

Ng i．oooo：oE oo 一一7．is；614E oo

≧1：量自9雛89．｝：孚面向善一3呈 3 1．gs7962E oo o．oeooeoE．ei HANPUKU KAIsU ＝ 4

X ＝， 1．4ST962E OO F（xl ＝ e．ooooooE−ol

Fig．6 Program and result of Newton−Raphson s method．

10−6とする。例題では，a＝20およびb＝10≧して4回の反復によって，x＝1．487962なる根が10一6以上の精度で求ま

っている。ところで，学習者の求めた根は，一〇．73〜1．9 の間に分布した。たいていの場合，反復回数は10回以下となっている。

方程式の解法としては，他に，von Mises法，はさみう

ち法，二分法，Homer法などを説明する。これらの解法

を全て説明した後，反復回数，誤差，プログラムの比較などを行う。

3．2連立一次方程式の解法

連立一次方程式には，構造材の応力解析，機械系の振動など工学上のいろいろな問題で定式化される。連立一次方程式の解法には，Gauss法， Gauss−Jordan法， Jacobi法，

Gauss−Seidel法， Dooli亡tle法， Crout法，．逐次加速緩和

法，共役勾配法など多くの方法がある。中でもGauss−

Jordan法は， Gauss法を改良した消去法と呼ばれる解法である。これは，計算原理が単純で枢軸要素の選び方を工夫することによって，相当精度よく求まる。すなわち，枢軸要素をce iステップでee i列の要素の中で絶対値最大のも

のを選択するように解を求めればよい。このような理由から，連立一次方程式の解法のプログラミングには，Gauss−

Jordan法を取り上げる。

n元連立一次方程式の係数と定数項から行列表示し，増

大行列Aを作る。その枢軸要素をaleleとすると，次の2つ

の手順を繰り返し，n回の枢軸操作によって得られた増大

行列の第〃＋1列の要素の値が，各変数のとる値となる。

（1）行列Aの第le行の各要素をalekで割る＠好＝α研α㈱）。

②行列Aの第i行の各要素からee le行の各要素1（ aileを乗じたものをそれぞれ引く（〆ゴゴ＝鰯一aika ki）。

この手順に従って，Fig．7のようなプログラムが作成できる。第2〜21行は，増大行列の読み込みおよび印刷をする

部分である。第22〜62行は，絶対値最大の枢軸要素の探

索，行の交換，掃き出し操作，その後の増大行列の印刷な

どを行なう部分である。第63〜68行は，求めた答を印刷する部分である。

求める連立方程式は次のようなものである。

｛

ax一 y 十 2 ＝＝4

x 十by十2a＝9

−X−2Y十。2 ＝＝2

ここで，係数の中のa，b， Cは次のように与える。すなわち，a＝（学籍番号の下2桁），う＝（出席番号），および。＝

（a＋b）。その他の係数は例題と同じ値を用いる。例題では，a＝4， b＝6， c＝5としてx・＝1， y・・r1，2＝1なる解が 10−6程度の精度で求まっている。学習者の求めた解は，x は0．012〜3．1の間に，yは0．11〜2の間に，2は0．019〜

6の間に分布した。

連立一次方程式の解法としては，他にGauss法， Gauss一一

Seidel法などを説明する。これらの解法を全て説明した

後，精度，演算回数，プログラムの比較などを行う。

3．3 最小二乗法

1組の変量（x，y）についてn個の観測値（Xi，） i），．（伝1，

2，……，n）が得られているとき， Xとyとの聞になんらか

の関数関係y＝f（X）があるとして，このデータにより良

く適合する曲線を求めたいことがよくある。このような場合，この関数f（x）を求めるのに最小二乗法を適用できる。

最小二乗法は，観測値と関数関係との差，すなわち近似のはずれの大きさの和を最小にすることによって実験式を得

ようとするものである。xとyとの関係は次のように表わ

される。

Jri ＝f（xi）＋Ei

ここで，εiは理論値との誤差である。この二乗の総和Σ

i−1

εi2を最小にするような関数ノ（x）を求める方法が最：小二乗

法である。この関数として，n次関数べき関数指数関

数，などが用いられている。べき関数および指数関数は変

換によって，一次関数に帰着できる。また，n次関数は，

(5)

数値解析における授業方法改善の試み宮地

L23456789012345678901234567890三Z34567890123今567890．1Z3↑567890LZ3今567890 1111111111222222222233533333334ら444ら44445555555555．66666666667

C 督砦骨 6AUSSoJORDAN 飼F＝THOD VIMENSION AtlO，11｝

blRITEC4，6000J

6000 FORMAT t lH ，，lt− GAUsS．・JORDANI ， t 養 OHou 菅管t曾ノノ事 C ．．． SHnKI DATA NO YoM：KoMl PEADC3，5000） N sqoo FORMATtl）r）

・N1＝fu＋1 00 20 r＝！，N

READc3，5Q2e）｛Atl，J），J＝1，Nl｝

5e2n FoRbelAT（7Flo．o｝

20 CONTrNUE

C ．．． SHOKI DATA Nn rNSATSU ptRITEtg，6020｝

6n20 FORMATcIH ，，・ KErsU GYORETsV，，

砦脚 A 1，J｝曾｝

oo 4e r＝t，N

urRITEt4e6040｝ tACT J），J＝1，Nl）

6∩40 FOf〜MA了ユ目，1P工lf】Oo2，

40 COAITIINUE c ．．． CnAuSs．JoRDAN HnLJ DV 30e 1・＝1，N TF｛1．EQ．N｝ GO TO lsO P昌。●o

ぐ o．● PIVOT NO SFNTA祥L／

DO 120 K＝1，N

IFcp．GE．．alt stA（k，ll）i・Gn To t？O P＝ABSCA（K．，1｝）

KP＝K

12n CONTINUE

rF（P．GT．1．eE 5t Ga TO 140 w RITE（4，6100）

610∩ FOR．MAT｛IHO，，x＋t IくONO MATR！X

c ．．．

140

00

6811

6120

c 瞬●，

2Zv

00

46

22

6220 280 300 c ．．．

菅

WAl，

1 HVND D…三SU 砦畳誉零｝

STOP

GYO NO KOKAN OO 16a J＝11Nt

i 1＝A｛1tJ）

A（1，J）＝AtKP，J）

A（KP，J）＝！一．，

CONTINUE

PIVOT；A（1，1）

WRITEt4，6120） r，T，PIVOT FQRIrA7「IIHO？脅PエVOT A 亀，12，9，弓■

12st）＝，tlPE15．6｝

HAKIDASHI SOUSA

DO 22◎ 」＝ヨ，卜しl A11，J，冨A ユ，」，ノPIVO丁

CONTINUE

Do 26e K＝1，pt1．

IFtK．EQ．r） GO To 260 AKI＝A｛K，1）

Do 240 J＝rfpq

ACK，J｝＝ACK，J｝．AK1 Alr，J｝

CONT．INUE

CONTrNUE

｛）o z80 敦；！，卜蓬

wRrTEt4，6220｝（AtK，J），J＝1，Nl）

FORMATCIH ，iPllEIO．Z｝

CONT！NUE CONTrNUE KOTAE NO INSATSU ．WR1TEt4，6300）

6300 FOR卜4A了｛！！IH ，10Xe璽畳 KOTAε 畳．，

wRITEc4，6320） utA｛r，Nl｝，1＝ltN｝

6320 FORMAT（IH IIOXt，Xct一：Z ）＝，，

｝ IPE15．6｝

STOP END

菅蒼蒼 GAU55繭JORPAN HOV 膏軸涛

KESSU

3．00E 4．eoE 2．ooE

GYO只ETSU OO 4gOO乞 oo s．onE OO 冑3・OOE PIVOT A（

1．oeE oo o・oo医噌01 0●00E−O乳

．b 1｝富 li25E

A｛r，J・）

oo ．9．ooE oo ．1．80E ol oo ．1．oaE oo ．s．oeE on OO 呂．00E OO 2の70E Oi

2．50E1el り5・50E OO

．onooooE oO Oe il．75E OO・．1．25E ＿3・75巨 00 曝1・ら3E 1・ユらE O真 2・95E p1VOτ へ監 Z， 2、雛瞬5騎5りOOO．0∈ OO

1．00E． OO O．OOE．．Ot o．oaE一一〇t

O．OOE．Ol 8．64E．Ol S．45E ユ。．oOE OO 冑2●09E OO 口5會36E o，ooE−el ．3．23E Oa ．1．29E

︷U11

000 OO1000

一次関数の拡張として考えることができる6このように，

最小二乗法による一次関数近似は最も簡単なものである

が，他の関数の基本となるものであるので，ここでは，この一次関数による実験式を得るζとを考える。関数ノ②＝

a＋bxとすると，正規：方程式を解くとaおよびbは次の式

から得られる。

123ら56789σ上Z3ら567890−234567890工23ら567890ユ23ら567890乱234︐．6789

11工11111．112222222222333ヨ333333ら4444444ら4555．5555555 01閉23456N

⊂ 黄黄鱒 LEAsT SOUARE，5 MFTHOD

（：玉」I KA帖SU K！酋JI

PIF，，ENS！ON X（IOO｝，YtlOO），FtlOD）

c 5の■ 」）A丁角 NO YOMr起。幼1 READ（3，5000｝ N SOOO FORMATCI5｝

READ（3，SelO｝（XCT），Y（1），1＝1，tsO

soio FORMATt6Flo．n｝

一RTTE（4，6000）

6σ0σ FORMAT 1臼，9X，8蒼斎膏正噸」： 5HrKτ e●

蒔．・klO A了EHAME 畳曽砦ノ！

一 IH II Nottgxtlxttlgxt．y−t 一 zax，tsurTEicH！e）

C ．．． 5EiKI HOUTEI．SH：KI NO KE！SU KE！SAN sUMX菖0・O

SUMY＝Q．O SU卜丹X2嵩O・O SUMY2＝O．O SUMXY＝O．O Do 100 1冨1，糊 SUMX＝5UMX＋XC！）

SUMY＝SUMY＋Ytl）

SUMX2＝SUMX2＋X（1）scve2 5UMY2胃SU凹Y2φY｛D巽聾2 5UMXY罵SU夙X．YゆX｛1｝営Y｛工、

！oo coNTrNuE

C ．．． KINJ！ 5HIKI NO KEISutA，B） NO KEISAN AN冨N

AVEXpSUMX／AN AVEY＝5VMYIAN

B罫 5Uト1XY鱈AN畳AVEX曇AVEY｝！

脅 15UMX2騨AN鱒AV〒…X管砦〜｝

A＝AVEY．AVExXe

RHO3 5UMXY口AN曇AVEX督AVEY｝！SORT ｛燭 SUMX2−AN鋳AVEX賛う」2、砦重5Uト4Y2 AN曇斎 AVEY畳費2｝｝

C ．一． 5UrTEr−CHr Te KETTEI KEISU NO KEisAN GI＝O．O

G2富O・O DO 210 1冨19卜．l F｛：，葺AゆB砦X｛：｝

Gl＝G1＋｛Ytr｝．FtT｝｝一一Z G2＝G2■lY｛1，甲AVFY，引引2 2iO cONTINuE

RR雛1．O噂（Sl／G2

C ．．． 5HOR1 KEKKA NO INsATSU

WR！TE（4，6320 ） tl，X（T），Ytl），Ftl｝，

韓．！魯t，N，

6320 FORMAT（IH ，r5，1P3E15．6｝

WRITE（4v6340） RHe，A，ReRR 6340 FORMAT（／／IHO，一 50UKAN−KEISU 一，／

費砦蒼弧帰

順畳膏

STOP

END

IH ，・ RHO 置曾g1PE15・6ノノ！

1日，璽曇 KAIKI HOUTE！旧9・

，5HrK！曇壁！

IH t Sf ＝ VtEIS．6t ，，＋ 99ε15の6t●， X■ノ！

IHO 黶CKETTEI KFIsU ＝，一EIS．6｝

leJl SHIscr ，NO ATEHAME

x

2．000000E 4．oooeooE 6−ooooonE s．ogoeooE 1．00000nE l．200000E．

111122000000

^Yr．．700000E−e1

3．4SOOOOE−O1 5．250000E−O1

7●099999E甲（》真 g．esooooE．．o1 1．105000E OO

畳砦曇

5UrTErCHI 1．595234E．O1 3．463806E．Ol S．33237SE．Ot

？．200950E−Ol 9．069S21E．Ot l．093SOgE OO

P！VO下 A 39 3｝冨一ヨ・22727ヨE OO 1．oeE oo

o．ooE．ol OoOOE飼01

O．OOE，一〇1 O．OOE−Ol 2．00E leoeE on o．ooE．al 3．oeE O．OoE1Ol 1．oOE Oo 4．OoE

000 000

斡 SOUKAN印KEIsU 曇 RHO ＝ 9 ．99e642E−el

蒼κoγ《ε

XlP 冨

xく2，鴇 x【3｝冨

曇

2．oooeotE oo 3．QooaooE eo．

4．ooooooE oO

Fig．7 Program and result of Gauss−Jordan s method．

曇 KA：．K峯 HOVT．E！甲5HIKI 輪 Y ＝ t ．2．733386E．02）． t

KETTET KEIsU i． 9．99S287E．O1

9．34286DE．03） X

Fig．8 Program and result of least square method．

(6)

＝xiyi−nxy

∂＝とゴユ

Zxi2−nT2

t＝tL

a＝アーう万．

ここで，アおよびアは，観測データのXおよびyに関する

平均値を表わす。最小二乗法で求めた実験式の実際の観測データにどれだけ適合しているかを知るために決定係数R を用いる。この一次関数の実験式を求めるプログラムは，

Fig．8のように作成できる。第4〜13行は，観測データを

入力し，表題を印刷する部分である。第14〜26行は，観測データの和などを求める部分である。第27〜36行は，係数

aおよびbを求め，相関係数を求める部分である。第37〜

45行ぽ1決定係数と推定値を求める部分である。第46〜57

行は，観測データ，推定値，係数α，bなどの印刷をする部分である。

実験式は，Table 2のデータに基づいて求める。各学習

者は，xの値を学籍番号の下2桁ずつ加えたものを用いる

ようにした。例題ではa＝一〇．0273，b＝0．00934， R＝0．999 なる値が求まっている。例題のデータも各学習者のデータもXが等間隔に変化するので，決定係数は同じになる。学習者の求めた係数αは，一〇．036〜0。023の聞に，bは0．023

〜0．18の間に分布している。

Table 2 lnput data for least square method．

20 40 60 80 100 120

O． 170 0． 345 0． 525 0． 710 0． 905 1． 105

最小二乗法では，この他に，二次関数べき関数，指数

関数などの適用について説明する。同じ観測データをこれらの関数の実験式をあてはめた場合どのように決定係数が異なるかなどを説明する。

3．4補間法

対数表，三角関数表などのように，xのn＋1個の点（Xe，…

…，κn）における関数値f（X）である観測値から，任意の独立変数に対する関数値を求める方法を補間法という。補間法には，1次補間法，2次補間法，Newton法， Lagrange法，

Aitken法， Aitken−Neville法など多くの方法がある。補間法の特殊な場合として，逆補間法および外挿法がある。中

でもLagrange法は，観測点（きざみ幅）が一定でも一定でなくても使用できること，観測点Xは順序通り並ぶ必要

がないこと，などの理由から他の方法より一般的である。

このような点から補間法のプログラミングには，Lagrange

補聞法を取り上げる。

n次のLagrange補聞多項式は，次のように表わされる。

P（x）＝ izn

鼈黶

Z（k40（x−xfe）／，go（xi−xk））ノω κ≒包 iCAft

この式から得られた値が，求めたい関数値の近似値とな

る。これをプログラム化したのが，Fig．9である。第2〜

14行は，観測データの読み込みおよび印刷をする部分である。第15〜21行は，任意の補間したい値の読み込みをする部分である。第22〜31行は，Lagrange補間多項式を分解して順次掛け合わせて近似値を計算している部分である。第 32〜34行は，補間値を印刷する部分である。

二2345678901234567890L234567890ユ23456

1ーユー正1111122222222223333333 ⊂ 畳曽鱒！N〒ERPOLAT＝ON 臼Y しAGRANGE ME7HOD

DIMENS！ON XtSO），Y（SO）

C 開● 5HOKI DATA NO YO隅エKOM＝

READC3，5000） N

5000 FORMAT（IS）

READt3vSO20） tXt！），Y（1），rtl，N｝

SO20 FORMATC6FZO．O）

C ．．． SHOK： DATA NO rNSATsU

WRITEC4，6000｝

600e FoRMATtIHOetac一一 LAGRANGE HOU Nr t

蒼，YORU HOKA陽畳畳砦●ノ！！

sc IH tsX NOetgX一一X e14X，tYe）

WRrTE（4，6020）（：，Xtl｝，y｛r），T＝1，NJ 6020 FoRMAT｛IH t2xtr5，IP2Els．6｝

c ．．． HoKAN suRu ATAr No yoM：KoMr

WRrTEt4，61001

6100 FORMAT ！ノ1H ，10X・19し HOKAN甲（：H1 曇・ノ聾 IH g16×9璽X●914X・身Y・，

loe READc3，sloe｝ xx SIOO FORMATtFIO．O）

rFcxx．EQ．o．el sTop C ．．． HoKAN−CHI NO KErSAN

p＝o．o

DO ！40 ！＝ltN

c＝yt！｝

Do 12e J＝1，N

：Ftl．EQ．J） GO To 120 c＝c一（xx．xtJ｝｝1（x（r）．xcJ｝）

ユ20 ⊂ONT工NUE ρ胃ρや（：

lgO CONTrNUE

（： ●●． HOKAN一（二H工 NO INSATsU

WRITEC4，6t20） XX，P 6120 FORMATCIH ，7X，IP2E！S．61 GO To 100

FND

賛管蒼 LAGRANGE HOU N塞 YORU HOKAN 管畳管

0乳7﹁∋456

N

…：書3898需匡 88 歪：亨88889蓬 89 ユ・8000QOE on 1・9000QOE O◎

曇 HOKAN，⊂HI 砦

x

Z．4SOOOOE OO 1．640000E oe 1．760000E oe

Y

2． 900eE eo 2．sssoooE no 2．756000E OO 2．924000E OO 3．egooooE no 3．2s6000E eo

ソ曇：自警1；呈1ε98 3・023883E ∩0 Fig．9 Program and result of Lagrange s method．

観測データは，Table 3の値を用いる。学習者は，最小

二乗法と同様にXの値を学籍番号の下2桁ずつ加えたも

のを用いるようにした。例題では，補間値は，x＝ 1． 45，

16．4，1．76について求め，それぞれ2．504，2．823，3．024

なる値が求まっている。

(7)

数値解析における授業方法改善の試み宮地

Table 3 lnput data for Lagrange s interpolation．

456789

¶←−←1111

2． 419

2．588

2． 759 2． 924

3．090

3． 256

補間法については，1次補間法，Newton法， Aitken−

Neville法などについても触れる。同じ観測データでどのような精度で求められるかを比較する。

3．5 数値積分法

実験値より積分を求めたいこと，関数が与えられていてもその不定積分が簡単に求められないこと，などは工学上の問題を取り扱う際によく遭遇する。このような場合数値積分が利用できる。数値積分には，台形公式，Simpson 1／3 公式，Simpson 3／8公式， Miine公式， Weddle公式， Hardy

公式など多くの公式があり，またその他にGauss積分，

Romberg積分などがある。中でも台形公式は実際の面積を

台形の面積で置き換えるというように原理が最も簡単で，

他の方法の基本になるものである。また，きざみ幅を小さくすることによって他の方法と同程度の精度を得ることも可能である。更に，偶関数，周期関数について言えば台形公式が最も良い積分法であることも知られている7）。このような理由から，数値積分のプログラミングには台形公式を取り上げる。

台形公式は，与えられた曲線y＝／（X）の面積を次の式で近似していく積分法である。

Si＝h（ノ（Xi）＋ノ（Xi＋1））／2

ここでhはきざみ幅である。従って，全区間の面積Sは，

区間の数だけこれらを加えればよい，すなわち，

s＝xsi

許1

によって求められる。この台形公式による数値積分のプログラムが，Fig．10である。第2〜14行は，積分区間の読み込み，そのデータおよび表題の印刷をする部分である。第 15〜22行は，．きざみの数を読み，積分値を求め，誤差計算して，それらを印罵する部分である。第25〜36行は，台形公式のサブルーチン副プログラムである。第37〜41行は，

積分する関数を定義する部分である。第42〜47行は，不定積分を定義する部分である。

積分する関数はア（x）＝exp（ex／のである。この関数において，a＝（学籍番号の下2桁）およびb・一（出席番号）として各学習者ごとに異なった関数の積分値を求める。積分区間は，〔0，1〕である。等分数は，6，12，36，60，120，498

1234567890123ら567890工2345678901234567890ユ234567

1111ニー1111222222Z2Z233﹂∋33333344444ら4 C 一一 TNTEGRAL BY TRAPEZOIDAL FORMULA

WRITEC4，6000｝

6000 FOR卜6AT tH gl膏黄砦 DAIKE！口50KU ．●

舜 INI YORU SEKIBUN 砦鹸蒼・，

10 READ｛3，5000） A，B soeo FoRMATt2Flo．o）

IFCA．EQ．B｝ STOP T￥eTF｛A，B）

WRITE｛4，60ZO） A，A，TY

6020 FORMATt／11H ，tSEKIeUN KUKAN （t，

sc IPEIO．3tl t−IEIo．3ft ）1／

蒼 1H ，，SH！N NO ATAI 30，F三5・6！ノ

ce iH ，t TOBUNSUtv4X，15EKrBl一

一 tUNcHr，一7xteGosA，1

C ．．． KAKV TOBUNSU N！ YOPU 5EKIBUN 10a READ（3，5tOO｝ N

5100 FOR団A了 15｝

！F（N．Eo．o） Go To ！e CALL DAIKEItA，B，N，S｝

E＝ABsCtTY．一S；／TY）

WRITE（4e6100） N，S，E 610e FORMATtlH ，3X，r5，IP3E15．6）

GO TO ユOO

ENP

⊂ 費費蒼 TRAPEZO！DA匙＿ FOR卜弓ULA

sUBROUTINE DAIKFItA，S，N，S）

H呂｛日鴨A，！FLOAT｛Nl

s＝o．o

Nl＝N．1 DO 20 1＝1，Nl

X冒A＋H畳FLOAT 1｝

5冨5やFξX，

20 CONTINUE

S＝H曇〔F（A⊃←F R｝＋2●O砦S，ノ2・O

RETURN ENDC XH KANSU NO TEIGr

FUNCTION FCX｝

F＝FxP｛23．0−X／9．0）

PETURN

ENF

C N管t SEκ1臼UN一⊂HI FVNCTION TF（A，B｝

TF； EXP 23．O管B／9，0｝齢EXP 23●0費A 糞ノ9．O，〕蒼9・O！23・O

RETURN

ENO

勇黄軸 DA！KE≡1層50KU NI YORU SEKIBUN 弥費曇

SEKI8VN KUKAN t O．OOOE−Ol ， 1．OOeE nO 1 5HIN NO ATAT ＝ 4．64AOfi8E OO

TOeUNsU 5EKIBUNCHI GOsA

追含：磯雛器 1：琴難1∈：9；

詣：：器灘881：塀留；1＝8：

120 4．64s252E oo 3．s2go21E−05 49B 4．64SO37E OO 1．Oe743iE−05

Fig．10 Program and result of trapezoidal formula．

である。例題では，a＝23およびb＝9としている。このよ

うな関数では等分数を増加すれば誤差は10−2〜10−5というように次第に減少している。学習者の求めた値もほぼ同様な傾向を示している。

数値積分としては，他にSimpson 1／3公式， Simpo且3／8

公式などにふれる。これらは，積分する曲線を直線，2次曲線，3次曲線で近似するものであるが，通常この順に精

度がよいといわれている。しかし，偶関数および周期関数などでは逆順に精度がよくなることを示して，数値解析の必要性を認識させるようにしている。

3． 6 常微分方程式の数値解法

電気工学，機械工学などの諸問題は常微分方程式で書き

表わされることが多い。それらの常微分方程式の一般解を

求あることは容易なことではない。また，現実の問題で

は，特殊解が必要であるので，数値解法が重要となる。常

微分方程式の数値解法は，Euler法，修正Euler法， Euler一

(8)

Richardson法， Euler−Gauss法， Runge−Kutta法， Runge−

Kutta−Gill法， Euler−Gauss予測子器正子法， Milne予測子

修正子法，Adams法， Hamming法など多くの方法がある。前進形解法で，最も単純で，誤差がきざみ幅の2乗程度であるEuler法について，フ。ログラミングの練習をす

る。

常微分方程式dy／dx一＝f（x， y）と初期条件（Xo， yo）が与えられて，xに対する関数値yを求める問題である。 Euler 法は関数値の初期条件を知って次の漸化式によって求める

ものである。

1234567890ユ234567890123ら567890三2345678901Z3ら567890ユ2

．111111111122222222223333333333邸舟南晦侮輪砧晦噂6555 C scNee nRDINARY DIFFERENTIAL EOuATrON

C EULER，S METHnD

WRITEC4，600e）

6000 FORMATCIH ，，一 t EULFR HOU，一砦・管砦涛・〃⊃

C ．．． SHOKT DATA NO YOMIKOMT OYOAr INSATSU READc3tsooo） xo，x1 tye

sooo FoRMAT（3Flo．e）

回R！丁E 4，6020， X（，，XP』gXO，YO 6n20 FORM，ATCIH SEKIRUN （VKAN （，t

sc IPEIO．3t， sEIO．31t ；tt t IH tlSHOKI−CH！ XO ＝tl 一 EIO．3，t NO TOK： t

ve 一Yo ＝ttFIO．3）_C

煙ｾ，EAD，3，，。、嚇．

SOIO FORMATt15）

！FCN．EQ．O｝ 5TOP H冨 XN−XO｝！FLOAT｛邸，

WRITEt4，60401 H

6040 FORMAT ！！1H gtKI2A団1 ：．，IPEIO・3ノノ

t IH tgXttXtt14XteYtf8Xt−5H！el

t tN No ATArtt7x GosA，1C ．．． SEKrSUN KErsAN x＝xo

Y＝YO Nl＝N＋1 DO 120 1＝1，Nl

TY＝TFCX）

E＝8AB5 TY印V，ノTYI WRITE−C4t6120） x，Y，TY，E

6120 FORMAT（IH ，IP4E：5．61 CALL EULER tX，Y，H）

X冒XO6FLOAT ：，9卜H

l20 CONTINVE GO了0 ユO

END

⊂ 聾費餐 EULER HOU

SU日ROU了rNE EUしER XgY，H，

Y＝￥＋H−FCXtY｝

RETURN

E閥D

C 膏費骨卜il騨SEK1BUN KANsU

FUNCTrON F（X，Y）

F＝t7．0−3．0−Y）t2S．O

RETURN END

C 鱒9卜曹 FUTE！ 5EK I8UN

FUNCTION TF｛X）

TF講7．O／300骨 1・0爾EXP 儒360！25・6tx3，

RETURN END

曇管麓 EULεR HOU 曽砦曇

§昌旨霞号ビ砦H誉UK日置 ≦ 8：888讐＝8葦，Nδ●孚82｝ oeo，冨 o・oooε＿Ol

K！ZAMr ＝

x

O．OOOOOOE−O1 1．250000E．O1 2．500000E．一〇1 3．7SOOOOE．Ol S．OOOOOOE．O1 6．2SooooE．el T．soooeoE．ol S．7500eoE−o1 1．000000E OO

Figユ1

1・250E胴0ユ

Y

O．OOOOOOE．O1 3．SOOOOOE 02 6．94750eE．02 1．034329E O1 1．368814E．O1 1．698281EeO1 2．022eO7E．O1 2．34246SE．Ol

Z・6517328∈騨〇二

5H1閥 NO A了Al O．OOOOOoE−O1 3．4T3850E−02 6．895992E−02 1．OZ6723E O1 1．358S26E−O1 1．685984E．O1 2．008271E−O1 2．325760E．O1 2．638S23E．Ol

GOsA

e．ooooooE−o1 7．S2747eE 03 7．469173E−O1 7．40TS99E−03 7．3999S7E−03 7．293861E−03 7．237960E−03 7．le2512E−03 7．126968E−03 Program and result of Euler s method．

Yi＋1＝Yi＋（Xi＋1−Xi）f（Xi， Yi）

これは，曲線を長方形で近似して積分を行なっていることになる。誤差は（Xi＿1−Xi）2程度である。この解法をプログラム化したのが，Fig．11である。第3〜14行は，表題の印刷，初期条件および積分区間の読み込み，その印刷などをする部分である。第15〜23行は，等分数の読み込みおよびきざみ幅の印刷をする部分である。第24〜35行は，漸化式を用いた関数値を求める部分である。第38〜42行は，

Euler法のサブルーチン副プログラムである。第43〜47行

は，関数ア（x，y）を定義した部分である。第48〜52行は，

特殊解を定義したものである。

常微分方程式dy／dx＝（a−by）／25をEuler法で解く。この方程式において，a＝（学籍番号の下2桁）およびb ・＝（出席番号）として，各学習者ごとに異なった常微分方程式を解く。初期条件は，x ・＝Oのときy＝oであるとする，積分区間は〔0，1〕である。等分数は，8，32，128，512である。例題では，a・＝＝7およびb＝3としている。この場合，等分数を増加すれば，誤差は10−2〜10−4というように次第に減少している。学習者の求めた値もほぼ同様な傾向を示している。

常微分方程式の解法としては，他にRunge−Kutta法，

Euler−Gauss予測子修正子法， Milne予測子修正子法などに触れる。これらの方法は，あるきざみ幅までは，次第に精度が良くなる。しかし，ある値を過ぎると繰り返し回数が増加しまるめ誤差が蓄積されて次第に精度が悪くなる。

これらのことをいろいろな関数について比較などして説明

する。

3．7モンテカルロ法

乱数は，標本調査などで対象を任意に抽出する場合と

か，現象のモデルを作って現象の様子を調べるシミュレーションを行なう場合などに用いられる。その乱数の発生の仕方にはいろいろあるが，電子計算機を用いてモンテカルロ法を行なう場合は，算術乱数を用いることが多い。このモンテカルロ法は，解析的な式が得られない場合，実験的な手法が不可能な場合などの科学技術上の複雑な問題を解くために用いられる。モンテカルロ法の適用例としては，

重積分，微分方程式，積分方程式，整数計画法，シミュレーションなどがある。ここでは，興味ある問題として解析的に正確な値を求めることが困難であり，また，他の数値積分法の適用が困難である重積分について，プログラミングしてみる。

問題は，次の積分値Sを求めることである。

s＝／iJ g・・・…fgf（xl， x2，・・・…，x．）dxldx2・・一…dx．

モンテカルロ法で，この値を求める手順は次のようであ

(9)

数値解析における授業方法改善の試み宮地

る6）。

（1）区聞〔O，1〕上のm個の一様乱数π1，u2，…，Umを発生し，

f1 一f（Ul，Pt2，…，Um）を求める。次に，区間〔O，1〕上の別の

一一rt乱数U 1，U 2，…，U mを発生し，．f2＝f（W 1，U 2、…，U m）を求

める。YJ下同様にして， f3，f4，…，fnを求める。

②積分値Sの推定値が次のようにして求まる。

S＝＝ fi／n i＝1

このモンテカルロ法による数値積分のプログラムがFig．12 である。第3〜14行は，積分区間の読み込み，そのデータおよび表題の印刷をする部分である。第19〜32行は，乱数を発生させ，積分値を計算し，10i2（i＝1，…，10）ごとにそ

ユ2345678901234567890上23456了890互234う6→890123456789012

1−1111﹁11112222222222333ヨ353555晦ら喬硲らり峰喬喬南555 （誉軸砦 NリMER！（Aし INTεGRATION

⊂ BY MONTE ⊂ARしO MFTHOO READ｛3，5000｝ A，g

500Ci FORMAT｛2FIO．O｝

TY＝TFtA，BI

WRrTEC4，6000） A，PsTY

6000 FORMAT CIH ，，kX MONTE CARLO HOU，，

一， NT YORU SuUCHI SE〈lt，

曇 ORUN ftscsc．！！〆

一 iH ，，5EKIFUN KUKAN Ci，！P

黄 E！O・3gサ，o，E工0・3g．｝甲ノ費 1H o●阜H！N NO ATAI 2曾，F15の6！ノ

一 1H ，tKAISU，，4X，，SFKIBUNCHI，，

一 ex

C ．．． 5HOKI−5ETTEI

GosAt｝

RN＝153．O N：O s＝e．e C ．．． SEKTBUN KEISAN DO 60 1＝ltle

20 DD 40 J＝ltlO URNI＝URAMCRptLl URN2＝URANCFIN1

5＝S＋F（URr．11eLIRKI2｝

4a CONTINUE

N＝N＋10

！F〔N・NE●10畳：考管2｝ 60 70 20 EKB＝SIFLOAT（N）

ERR＝ABsCCTY−5KM｝11Y）

WRrTEt4｝6i（le｝ N，gKR，ERR 6100 FORMATtlH ，15，IP？F15．6｝

6e CONTINUE 5TOP FND

C 聾誉聾 GENEF〜ATION OF RA吋DOM NUM臼ERs FUN⊂T工ON URAN｛R卜』，

RN＝AMOD｛765．OscRN，3276S．O）

URAN＝RNI32768．o

RETURN END

C曇曇骨 H1SEKIE〕UN KA卜ISU NO TEIGI FVNCTION F（X，Y｝

F：X費ヒX曇管2ゆ2●0｝畳曇！06NY砦EXP Y／73●0｝

RETURN ．

END

C 鱒黄管 5H1N NO SEK！RUN⊂H：

FVNCTION TF（A，B）

TF霧｛ B誉膏2←2●0｝砦砦2●5閏 A畳藷2φ2．O｝賛曇2．5 管】ノ5・0皆 EXP B！73■0，砦 B疇7ヨ00｝o EXPtAIT3en）vetA．73．ol｝t73．O

RETURN

END

督畳骨 MONTE CARLO HOU NI YORU SUUCHI sEKIBUN 碁畳脅

SEKIeVN

SHIN NO KA！

Pと

紹

碧8

2gg

自皇8

！ooo

KUKAN く ATA：胃

SEKTBUNCHI

ユ．2ユ95ら3E OO 1．047086E oO 9．229667E．O1 8．699827E．O1 9．272032E−Ol 工の003631ε Oh 1．030202E oO 1．04S66？E OO 1．024394E oo i，003443E OO

O．OQOE−O1 1．oo22goF

， 1．000E OO ）

oo GOSA

2．167692E−O1 4．470408E．02 7．913285E−02 1．319963E．O1 7．490604E．02 1．34B282E．03 2．7S5888E 02 4．32B800E．02 2．206349E−02 1．160122E一一〇3

Fig．12 Program and result of Monte−Carlo method．

の結果を印刷する部分である。第35〜40行は，一様乱数を発生させるサブルーチン副プログラムである。第41〜45行は，．積分する関数を定義する部分である。第46〜52行は，

不定積分を定義する部分である。

モンテカルロ法で解く重積分は次のようなものである。

S一：轤煤轤№?i・・＋の・・…xp（・／の鋤

この重積分において，a＝（出席番号）およびb一＝（学籍番号

の下2桁）として，各学習者ごとに異なった常微分方程式

を解く。積分区間は〔0，1〕である。関数の計算回数は1000

回までとしている。例題では，θ躍2およびb＝73として

いる。真の積分値は，1．002280である。誤差は計算回数が増加するにつれて10−1〜10−3というように次第に減少している。学習者の求めた積分値は，0．931〜67．3の間にある。

誤差については，例題と同様な傾向を示している。

モンテカルU法としては，算術乱数の発生法など説明

し，他に待ち行列，在庫管理，Buffonの針などについて触

れる。

3．8 シミュレーション

シミュレーションは，前項3．7で述べたモンテカルm法

で行なうものもある。しかし，ある物理モデルを設定し，

そのモデルの運動を表わす微分方程式を立てて，モデルの定数をパラメータとする。このパラメータを変えることによってそのモデルがどのようにふるまうかを見るようなことも現実問題としては多く見られる。ここでは，このよう

なシミュレーションを考える。その1つの例として機械工

学でなじみ深い振動系の基礎的なものとして減衰自由振動を取り上げる。

1自由度機械振動系の減衰自由振動の状態方程式は次の

ように表わされる。

dx ．＝v dt

窪一÷瀞一・・9・（・）蕩

ここで，x， v， t， m，ゐ， C， fはそれぞれ変位，速度，時

間，質量，バネの弾性係数，粘性減衰係数，摩擦力であ

数値解析における授業方法改善の試み

数値解析における授業方法改善の試み

宮 地 潜＊

（昭和54年5月2日受理）

Improvement of Teaching Method in Numerical Analysis

1．序

ここでは，著者の担当する数値解析の授業を工夫し，一応 の成果をあげていることを報告し，情報処理教育および教 育方法改善の参考に供したい。

＊機械工学科

プログラミングの好適題材として数値解析がある。しか し，本校では従来数値解析はプログラミングに余り利用さ

ングおよびそれぞれの問題に合った計算方法に習熟するよ うに考えた。すなわち，数値解析の各種の手法を説明する

し，プログラムを作成させる演習を行なうように授業を展

者の固有の学籍番号，および出席番号などを用いて変えて

に比し積極的に演習する姿が多くみられ，効果がかなりあ るように思われる。

以下では，まず授業方法の概略を説明する。次に，数値

解析の授業で取り上げている内容，．すなわち，方程式の解

法，連立一次方程式の解法，最小二乗法，補間法，数値積

分法，常微分方程式，およびモンテカルロ法の代表的解法

を1．つずつ取り上げどのように演習を行なっているかを順

Solution

equation

Simultaneous linear equation

Least square

Interpolation

Numerical integration

differential

equation

Monte Carlo

Contents

Beirstow s method

von Mises method

Bisection method

M・・…P・・at・・n｛1鵠尊詠Mpy）

inverse matrix

Gauss−Seidel s method

翻黙黙翫）

Lagrange s method Newton s method Aitken s method Aitken−Neville s method

Gauss method

Modified Euler s method Euler−Gauss method

Runge−Kutta−Gill s method

＊Milne s method

11⊥11 444一4

2

Steps

49 52 134 52 T9 T5 T4 T5 T6 U9 ﾐ

117

Data linage

3229臼ウ臼

4

No． of presenters

75一75

69 6FO4 K＝一46扉

62

Rate of completion

87︸84

90 W4 X6 W9 W2

器＝﹁93︸96

FひFOO

jobs

32

弓＝至一β 22． 1

2

次説明する。最後に，演習結果のまとめを示し考察を加え．

2．授業方法の概説

カード．ディックを編成することにしている。ただし，Fig．

3のようなパンチカードで編成を希望する者には，．演習の 日までにカードをパンチしておくことを条件に許可してい る。コントm一ルカードなどを付けて，カードディックを

Modification

no

yes Producti6n

yes Grading Stop

analysis．

数値解析における授業方法改善の試み 宮 地

Fig．2 Mark card for preparing FORTRAN programs．

● 鵬 Fig．3 Punch card for preparing FORTRAN programs．

プログラムを修正すると同時にLP用紙とチェックカード を提出し，再び計算依頼することになる。

課鵬号感激 番雛 回

考 察

また，数値解法の区切りごとに，精度，収束回数，など について各解法を比較し，解説している。そして，数値計 算における注意すべきことを説明するように努めている。

数値解析の授業で扱っている内容は，次の8つに大きく 分けることができる。（1）方程式の解法，｛2漣立一次方程式

分方程式の数値解法，（7）モンテカルロ法，⑧シミュレーシ

ョン。この節では，これらの内容の中から，それぞれ代表 的な解法を1つずつ取り上げて演習方法を少し詳しく説明

宮地 ^潜＊

ここでは，著者の担当する数値解析の授業を工夫し，一応の成果をあげていることを報告し，情報処理教育および教育方法改善の参考に供したい。

プログラミングの好適題材として数値解析がある。しかし，本校では従来数値解析はプログラミングに余り利用さ

ングおよびそれぞれの問題に合った計算方法に習熟するように考えた。すなわち，数値解析の各種の手法を説明する

に比し積極的に演習する姿が多くみられ，効果がかなりあるように思われる。

3のようなパンチカードで編成を希望する者には，．演習の日までにカードをパンチしておくことを条件に許可している。コントm一ルカードなどを付けて，カードディックを

数値解析における授業方法改善の試み宮地

プログラムを修正すると同時にLP用紙とチェックカードを提出し，再び計算依頼することになる。

課鵬号感激番雛回

考察

また，数値解法の区切りごとに，精度，収束回数，などについて各解法を比較し，解説している。そして，数値計算における注意すべきことを説明するように努めている。

数値解析の授業で扱っている内容は，次の8つに大きく分けることができる。（1）方程式の解法，｛2漣立一次方程式

ョン。この節では，これらの内容の中から，それぞれ代表的な解法を1つずつ取り上げて演習方法を少し詳しく説明

技術関係の問題を解くとき，方程式の根を求めることがよくある。方程式の解法には，Cardano法， Ferrari法，

New㎞一Raphscm法によるノ（X）＝Oの根を求める漸化式は，Taylor級数のxに関して1次までの項によって求めら

の印刷および初期値の読み込みなどの準備部分である。第 11〜21行が反復計算して根を求める部分で，1回計算することに得られた近似値を印刷するようにしてある。第22〜

27行は，求められた根を印刷する部分である。第30〜34行は，根を求める関数を定義した部分であり，第35〜39行は，その導関数を定義した部分である。

根は，超越関数f（x）＝・ex2−beXについて求める。この関数において，a＝（学籍番号の下2桁）およびδ＝（出席番号）

10−6とする。例題では，a＝20およびb＝10≧して4回の反復によって，x＝1．487962なる根が10一6以上の精度で求ま

っている。ところで，学習者の求めた根は，一〇．73〜1．9 の間に分布した。たいていの場合，反復回数は10回以下となっている。

を全て説明した後，反復回数，誤差，プログラムの比較などを行う。

連立一次方程式には，構造材の応力解析，機械系の振動など工学上のいろいろな問題で定式化される。連立一次方程式の解法には，Gauss法， Gauss−Jordan法， Jacobi法，

のを選択するように解を求めればよい。このような理由から，連立一次方程式の解法のプログラミングには，Gauss−

②行列Aの第i行の各要素からee le行の各要素1（ aileを乗じたものをそれぞれ引く（〆ゴゴ＝鰯一aika ki）。

この手順に従って，Fig．7のようなプログラムが作成できる。第2〜21行は，増大行列の読み込みおよび印刷をする

どを行なう部分である。第63〜68行は，求めた答を印刷する部分である。

ここで，係数の中のa，b， Cは次のように与える。すなわち，a＝（学籍番号の下2桁），う＝（出席番号），および。＝

く適合する曲線を求めたいことがよくある。このような場合，この関数f（x）を求めるのに最小二乗法を適用できる。

最小二乗法は，観測値と関数関係との差，すなわち近似のはずれの大きさの和を最小にすることによって実験式を得

数値解析における授業方法改善の試み宮地