今回のテーマ: 数値積分
前回から演習の後半戦で、物理現象の簡単な数値シミュレーションの入口にあたる実習を行っている。 その主旨は、 手計算でやるのが大変な問題を、コンピュータに代わりにやってもらう ことであるともいえる。 前回はその手始めとして、非線形方程式の数値解法を扱ったが、 今回は 関数の定積分の計算を、コンピュータに代わりにやっててもらうことを考える。 定積分の計算は、科学技術計算のいろいろな場面で登場する。 その上、多くの場合はそもそも 手計算で解くことが極めて困難 であるため、どうしてもコンピューターに頼らざるを得ない問題である。
台形則による数値積分 (その1)
例題として、関数 \( f(x) \) を \[ f(x)=\dfrac{1}{1+x^2} \] ととり、これを区間 \( 0 \le {x} \le 1 \) で積分した値 \[ S\equiv\int_0^1f(x)dx \] を、台形則を使って求めるプログラムを作成してみよう。 必要な手順 (の大部分) を以下に示す。 ただし例の如く、肝心なところは (わざと) 抜いている ので、各自で考えて埋めてほしい。
また以下の「ひな型」では、個々の「短冊」の面積を「台形の面積の公式」から求め、その総和をとる方法になっている。 もちろんこの「ひな型」通りの方法ではなく、講義で示した「台形則」の最終的な表式に則ったプログラミングでもOK (その「ひな型」は ここ から)。
なお、このプログラムのミソ (というか動作の流れ) は以下の通り。
-
定積分の値を計算するのに必要な「器」を宣言する。 例えば「総和」をとるのに使う souwa、「短冊」の幅を表わす dx、及び定積分の値を入れる sekibun を用意する。
「短冊」の数や、その両端での x や y の値、および「短冊」の面積を入れる「器」を宣言する。 ここでは、
「短冊」の数 を表わす整数型変数 kosuu を定数 (パラメータ) として宣言 し、同時にその値を設定 している。 parameter という「属性」の詳しい説明は ここ を参照のこと。
「短冊」の両端での x や y の値、および「短冊」の面積 を入れるための 実数型配列 を宣言し、同時に定数 kosuu を使って配列の大きさを指定 している。 なお、配列の添字の範囲の指定を (0:kosuu) とすると、その配列の要素は 0番めから kosuu 番め として参照され、要素の数は kosuu+1 個となる。 配列の宣言のしかたの詳しい説明は ここ を参照のこと。
整数型の定数を用いて配列の大きさを指定する ことにより、その 定数の値だけを変更するだけで、配列の大きさに関わる箇所全てを一度に変更 することができるという利点がある (もちろん、「そうなるようなプログラム」になっていれば、の話だが)。
定積分を計算する区間を指定するのに必要な「器」を宣言する。 ここでは「始点」x_kaishi と「終点」x_shuryo の双方とも定数 (パラメータ) として宣言し、同時にその値を設定している。
定積分の値を計算するのに必要な変数に値を入れる。 当然ながら、「短冊」の数が必要な時は定数 kosuu を使って指定 するようにすること。 同様に、積分区間に関する情報が必要な時には定数 x_kaishi と x_shuryo を使って指定 するようにする。
個々の「短冊」の面積を計算し、その値を実数型配列 s のしかるべき場所に入れてやる。 「台形則」というからには、「短冊」の形が台形だと思って面積を計算 してやればよい。
定積分の値を計算し、変数 sekibun に入れる。 これまでに (しつこいほど?) やった方法で、個々の「短冊」の面積の総和をとればOK。
計算結果を表示する。 ここでは結果の正しさを判断しやすくするため、sekibun そのものの値だけでなく、それを4倍した値も同時に表示するようにしてある。
| p | r | o | g | r | a | m | e | n | s | h | u | 1 | 2 | a | |||||||||||||||||||||||||||||||
| i | m | p | l | i | c | i | t | n | o | n | e | ||||||||||||||||||||||||||||||||||
| i | n | t | e | g | e | r | , | p | a | r | a | m | e | t | e | r | : | : | k | o | s | u | u | = | 1 | 6 | |||||||||||||||||||
| r | e | a | l | , | d | i | m | e | n | s | i | o | n | ( | 0 | : | k | o | s | u | u | ) | : | : | x | , | y | , | s | ||||||||||||||||
| r | e | a | l | , | p | a | r | a | m | e | t | e | r | : | : | x | _ | k | a | i | s | h | i | = | 0 | . | 0 | , | x | _ | s | h | u | r | y | o | = | 1 | . | 0 | |||||
| r | e | a | l | : | : | s | e | k | i | b | u | n | , | d | x | ||||||||||||||||||||||||||||||
| i | n | t | e | g | e | r | : | : | i | ||||||||||||||||||||||||||||||||||||
| d | x | = | ( | x | _ | s | h | u | r | y | o | - | x | _ | k | a | i | s | h | i | ) | / | r | e | a | l | ( | k | o | s | u | u | ) | ||||||||||||
| d | o | i | = | 0 | , | k | o | s | u | u | |||||||||||||||||||||||||||||||||||
| x | ( | i | ) | = | x | _ | k | a | i | s | h | i | + | d | x | * | r | e | a | l | ( | i | ) | ||||||||||||||||||||||
| y | ( | i | ) | = | 1 | . | 0 | / | ( | 1 | . | 0 | + | x | ( | i | ) | * | * | 2 | ) | ||||||||||||||||||||||||
| e | n | d | d | o | |||||||||||||||||||||||||||||||||||||||||
| d | o | i | = | 1 | , | k | o | s | u | u | |||||||||||||||||||||||||||||||||||
| s | ( | i | ) | = | |||||||||||||||||||||||||||||||||||||||||
| e | n | d | d | o | |||||||||||||||||||||||||||||||||||||||||
| s | e | k | i | b | u | n | = | 0 | . | 0 | |||||||||||||||||||||||||||||||||||
| d | o | i | = | 1 | , | k | o | s | u | u | |||||||||||||||||||||||||||||||||||
| s | e | k | i | b | u | n | = | s | e | k | i | b | u | n | + | s | ( | i | ) | ||||||||||||||||||||||||||
| e | n | d | d | o | |||||||||||||||||||||||||||||||||||||||||
| w | r | i | t | e | ( | * | , | * | ) | k | o | s | u | u | , | s | e | k | i | b | u | n | , | s | e | k | i | b | u | n | * | 4 | . | 0 | |||||||||||
| e | n | d | p | r | o | g | r | a | m | e | n | s | h | u | 1 | 2 | a |
これを基にして、台形則を使って上記の定積分の値を求めるプログラムを作成せよ [みほん]。 このプログラムを (例えば) enshu12a.f90 という名前で保存し、例の如く gfortran enshu12a.f90 とした後に、
$ ./a.outと実行し、期待通りの答 (もちろん \(\pi/4\) です) が得られているか確認せよ。 正しく動いていれば、だいたい以下のような出力が得られるはずである。
16 0.785235345 3.14094138
正しく動いていることが確認できたら、本日の演習の出席確認のため、このプログラムを 末尾のフォーム 経由で提出せよ。
台形則による数値積分 (その2)
上で作ったプログラムで得られた計算結果は、真の答 (\(\pi/4\)) にさほど近いものではなかったはずだ。 しかし当然予想されるように、「短冊」の幅を小さくすれば、言い変えれば「短冊」の数を多くすれば、真の答により近い答が出てくる はずである。 そこで、「短冊」の数と、計算値の正確さの間の関係 を確かめてみよう。
上で作ったプログラムのうち、計算結果の出力部分を以下のように変更してコンパイル・実行し、真の答からのずれの大きさを見積ってみよ。 なお、atan(1.0) とは、プログラム中で手っ取り早く \(\pi/4\) を与えるための仕掛け である。 (atan とは逆正接関数と呼ばれるもので、tan の逆関数である。ここでは \(\tan(\pi/4)=1\) であることを利用して、\(\pi\) の値を与えている。)
| w | r | i | t | e | ( | * | , | * | ) | k | o | s | u | u | , | s | e | k | i | b | u | n | , | s | e | k | i | b | u | n | * | 4 | . | 0 | , | a | b | s | ( | s | e | k | i | b | u | n | - | a | t | a | n | ( | 1 | . | 0 | ) | ) |
次に、「短冊」の数 kosuu を 2 から 512 あるいは 1024 程度まで2倍ずつ増やしながら同じ計算 を行い、計算結果の正確さがどれくらい増していくかをを調べてみよ。 具体的には、kosuu を倍に増やすと、真の答からのずれがどれくらい小さくなってくれるか を見積ってみよ。 さらに kosuu と真の答からのずれを両対数グラフにプロットしてみると、どのような傾向が見られるか [みほん(その1) または (その2)]。
台形則が正しく動いていれば、「短冊」の数を倍にする、あるいは「短冊」の幅を半分にすると、求まる答の誤差は約1/4に減少する ことが見てとれるだろう。
余力のある人向けのおまけ: 自前の「関数」を定義する
これまでのプログラミングでも、例えば「絶対値をとる」操作や「平方根をとる」操作を行うために、いくつかの「関数」を使ってきた。 これらは 組み込み関数 と呼ばれるもので、よく使われるが故に (Fortran に限らずプログラミング言語の中に) あらかじめ用意されているものであった (他の例は ここから)。
こうした組み込み関数の他にも、プログラマが自分用の関数を新しく定義 し、その自前の関数を使ってプログラムを書くこともできる。 ここ に挙げたプログラムは、自前で関数を「副プログラム」 (subprogram) として定義し、その関数を使って計算させる手順の例を示している。 詳しい説明は ここ を参照のこと。 ここでは Fortran 90 の「内部副プログラム」という形式で自前の関数を my_kansuu という名前で定義し、ある値 (副プログラムの中では z) を入力すればこの関数 my_kansuu がしかるべき値を返すように設計してある。 ちなみにこのプログラムは、第10回 めの課題の1つである \(x^2=2\) を二分法で解くプログラムであり、この中で定義した自前の関数 my_kansuu はこの式の「左辺」-「右辺」を計算するものである。
余力のある人は、今日の課題である台形則による数値積分のプログラムを、自前の関数を定義して利用するスタイルで書き直してみよ。 具体的には、被積分関数 \( f(x) \) を自前の関数として定義し、この関数を使って「短冊」の両端での y の値を与えてやればよい [みほん]。
また 第8回 で紹介した サブルーチン (subroutine) も、Fortran で用意されている「副プログラム」の一種である。 ここ に挙げたプログラムは、サブルーチンの使い方の例であり、上と同じ問題をニュートン法で解いている。 サブルーチンの詳しい説明は ここ を参照のこと。
シミュレーション業界に限らず「実用的な」プログラムは、たいがいが (何らかの形で) 副プログラムを利用して作成されている。 プログラムは多くの処理の組み合わせから構成されているのだが、複雑なプログラムでは必要な処理の数が膨大になってしまい、「整然とした」プログラムに仕上げるのはそう簡単にはいかない。 そこで、「単一のプログラムを書く」という発想を最初から放棄してしまうのが一般的である。 その代わりに 処理をいくつかの小さな単位に分割 し、そして それぞれの処理を別々にプログラミング しようという考え方をとる。 こうすることで、プログラム全体の構造がすっきりする だけでなく、プログラムの 変更・修正・テストがやりやすくなる という利点がある。
台形則による数値積分 (その3)
ここで作ったプログラムを利用して、以下の定積分を計算してみよ。
- 高校の微積分で出てきた (はずの) 定積分。 \( \displaystyle\int_{0}^{1}\sqrt{x}dx \) とか \( \displaystyle\int_{1}^{2}\frac{1}{x}dx \) とか。
- 平均値が50、標準偏差が10の正規分布を表わす関数
\( \displaystyle G(x)=\frac{1}{10\sqrt{2\pi}}\exp\left[-\frac{(x-50)^2}{200}\right] \) の定積分。
例えば
- \( \displaystyle\int_{40}^{50}G(x)dx \) とか \( \displaystyle\int_{50}^{60}G(x)dx \) とか
- \( \displaystyle\int_{30}^{50}G(x)dx \) とか \( \displaystyle\int_{50}^{70}G(x)dx \) とか
- 三角関数に関する定積分。
これらは手計算でもできるようになっておいてほしい。
特に2番め以降の「直交関係」は「フーリエ級数展開」で頻出する。
- \( \displaystyle\int_{0}^{\pi/2}\sin{x}dx \) とか \( \displaystyle\int_{0}^{\pi/2}\cos{x}dx \)
- \( \displaystyle\int_{0}^{2}\sin^2(\pi{x})dx \) とか \( \displaystyle\int_{0}^{2}\cos^2(\pi{x})dx \) とか
- \( \displaystyle\int_{0}^{2}\sin(\pi{x})\cos(\pi{x})dx \) とか \( \displaystyle\int_{0}^{2}\sin(\pi{x})\sin(2\pi{x})dx \) とか
余力がある人向けのおまけ: 中点則
今日はじめに紹介した台形則による数値積分では、積分区間を細かく分割して作った個々の「短冊」を「台形」だとみなし、その面積の総和をとることで定積分を求めてきた。 しかし数値積分の本質的なところは、「短冊」の面積の総和をとることにあるので、個々の「短冊」の面積は (適切でありさえすれば) どのように求めても OK のはずである。 ということで、「短冊」の面積の計算のしかたをとり替えることで、異なった数値積分法を導くことができる。 以下では台形則と並んで基本的 (だと亀山が思う) な数値積分法として、「中点則」(midpoint rule) を紹介する。
台形則の場合と同様に、\({x_{i-1}}\le{x}\le{x_{i}}\) を占める短冊の面積を \(S_i\) と書くことにする。 中点則では、個々の「短冊」の面積 \(S_i\) を以下の表式により求める。 \[ S_i=(x_{i}-x_{i-1})\times{f(x=\frac{x_{i-1}+x_{i}}{2})} \] 即ち、個々の「短冊」を 幅が \(h\) (\(=x_{i}-x_{i-1}\)) 、高さが \(f(x=\displaystyle\frac{x_{i-1}+x_{i}}{2})\) の長方形 だとみなして面積を求めていくのが中点則である。
余力のある人は、中点則に基づく数値積分のプログラムを作成・実行し、期待通りの答が得られていることを確認せよ [みほん]。 さらに余力のある人は、「短冊」の数 kosuu を 2 から 512 あるいは 1024 程度まで2倍ずつ増やしながら同じ計算 を行い、中点則による計算結果の正確さがどれくらい増していくかをを調べてみよ [みほん]。 この手順によって、「短冊」の数を増やしたことによって正確さの向上する速さ を見積もることができるのだが、その速さを中点則と台形則で比較してみよ。
(講義ではやってないけど) 進んだ話題: 「乱数」を使った定積分の計算
2回生後期の「地球物理学実験」で考えた「乱数を使ったゲーム」を応用することで、定積分の計算を代用することもできる。 もちろん通常はこのようなやり方をする必要はないのだが、「止むに止まれぬ場合もある」ということで。 なお、乱数を用いて数値計算やシミュレーションを行う方法は一般に 「モンテカルロ法」 と呼ばれ、この名前はカジノで有名なモナコ公国の地区の名称 Monte Carlo に由来している。
こうした定積分の計算の簡単な (かつ多くの教科書で扱われている) 例は、円周率 \(\pi\) の計算である。 2年生後期の「地球物理学実験」で (ひっそりと?) 紹介されたのを覚えている人もいるだろうが、その手順は以下の通りである。
- 2N 個の乱数 \(r_1,\cdots,r_{2N}\) を、0 から 1 の範囲で発生させる。
- これらの乱数を用いて \(x_i=r_{2i-1}\) かつ \(y_i=r_{2i}\) とし、これらを x座標、y座標とする N 個の点を考える。
- これらの N 個の点のうち、原点からの距離が1以下 の点の数を数える。 要するに、いま考えている点 \((x_i,y_i)\) が、原点を中心とし、半径1の4分円の中にあるかどうか? をチェックしている。
- 上の条件を満足する点の割合として、\(\pi/4\) の近似値が求まる。 要するに、この手順で 一辺が1の正方形の面積と、半径1の4分円の面積の比 を近似的に見積り、これによって \(\pi/4\) の近似値を求めている。
もちろんモンテカルロ法を使って「普通の」数値積分をすることだって可能である。 例えば、\(\displaystyle\int_a^bf(x)dx\) という定積分をモンテカルロ法で計算するには、以下のような手順をとればよい。
- N 個の乱数 \(x_1,\cdots,x_N\) を積分区間 \({a}\le{x}\le{b}\) の範囲内で発生させる。
- 各々の乱数 \(x_i\) の値における \(f(x_i)\) の値を計算する。 要するに、\({a}\le{x}\le{b}\) の区間での \(y=f(x)\) というグラフの形状を、\(y=f(x_i)\) という「真っ平ら」な関係で置き換え てしまおうとしている。
- 幅が \(b-a\) (積分区間の幅)、高さが \(f(x_i)\) の長方形の面積を考える。 要するに、いま考えている長方形の面積で、求めたい定積分の値を代用 させてしまおうとしている。
- 各々の乱数に対応してできる N 個の長方形の面積の平均 (というか期待値) として、定積分の近似値を求める。
余力のある人は、モンテカルロ法による数値積分の計算を実行してみよ。 さらに、発生させる乱数の数を増やしながら同じ計算 を行い、真の答からのずれがどう変化しているか を調べてみよ。 また十分 (例えば小数点以下4桁程度まで) 正しい答を求めるには、どの程度の数の乱数を用いてやる必要があるか。
この例でみるように、モンテカルロ法による数値積分には、発生させる乱数の数を増やしてもなかなか答の正確さが向上してくれない という欠点がある。 しかしその反面、何重の積分であっても全く同じやり方でOK という (ありがたい) 性質があるため、多重積分 をしたい場合の「最後の手段」として使われることがある。
出席確認用プログラム提出フォーム
Moodle のコースページ 以下の 課題ファイル提出 により、各自で作成した (enshu12a.f90 に相当する) 台形則を用いて \(\displaystyle\int_0^1\frac{1}{1+x^2}dx\) という定積分の値を計算する Fortran 90 プログラムの最終版を提出せよ。