２００３年度・数理解析・計算機数学２・第４回 1
● 今後の講義予定（変更）
第４回(11/04)
【Ｃ言語及びその他(1)】
•関数ポインタ
【浮動小数点演算(4)】
•ニュートン法の収束
•逐次近似 第５回(11/11)
【浮動小数点演算(5)】
•多角形近似によるπの計算と数値誤差
•加速 第６回(11/18)
【浮動小数点演算(6)】
•級数和の計算 第７回(11/25)
【常微分方程式の数値解法(1)】
•常微分方程式の数値解とは何か
•オイラー法とその誤差 第８回(12/02)
【常微分方程式の数値解法(2)】
•改良／修正オイラー法
•数値的不安定性 第９回(12/09)
【常微分方程式の数値解法(3)】
•ルンゲ・クッタ法
•ルンゲ・クッタ法の安定性
•予測子・修正子法 １２月１６日：休講（予定）
第１０回(01/06)（補講）
【行列計算(1)】
•行列の和・積・スカラー倍 第１１回(01/13)
【行列計算(2)】
• （詳細未定）
第１２回(01/20)
【行列計算(3)】
• （詳細未定）
第１３回(01/27)
【行列計算(4)】
• （詳細未定）
● 前回の講義のまとめ
【浮動小数点演算の誤差】
1.浮動小数点数の四則演算（加減乗除）には演算誤差が生じる可能性がある. 浮動小数点数a, b∈Fに対して,その演算誤差は
|(ab)−(a·b)| ≤2εM|a·b|
と評価できる.ここで, “·”は加減乗除のいずれかであり,a·bは,a,b∈Fに対する正しい演算 結果,abはその演算結果をFに丸めたものである.
すなわち,Fの四則演算は相対誤差2εMを持つ.
2.特に,「減算」に関しては,少なくとも１桁の「保護桁」が必要であり,保護桁なしの減算の相 対誤差は極めて大きくなってしまう.
3.初回の例an= (1 + 1/n)nの計算においては,xn= (1 + 1/n)をFで表現するための誤差εM
(|xn−xn| ≤(1 +εM)|xn|)から積算される誤差neδだけでなく, (xn)nを求めるための（例え ば）n−1回の乗算による誤差(n−1)eεMも計算結果に組み込まれる.したがって,誤差の総 和は|e−an| ∼C(1/n+CnεM)程度となる.
ex04.tex,v 1.3 2003-11-02 15:59:52+09 naito Exp
2 ２００３年度・数理解析・計算機数学２・第４回
【数値の計算】
1. 目的の数値αが含まれる区間Inを,順に２分割して区間を縮小して近似列を作る,「区間縮小 法」（「二分法」）.
2. 区間縮小法では,近似値xnを「区間の中央値」と定めると,
|In|= 2−n|I0|, εn=|xn−α|<1
2|In|= 2−(n+1)|I0|, が成り立つ.特に,絶対誤差δでαを求めるための繰り返し回数は
εn<2−(n+1)|I0|< δ で求めることができる.（α=√
2の時には,I0= [0,2]ととればよいので,n >log2δである.）
3. 区間縮小法は,「左右どちらの区間に真の値が入るかを判定する関数」f(x)に関する,中間値 の定理を適用した方法と考えることができる.
4. 滑らかな関数f の零点ξ(f(ξ) = 0)を求める時,f がξの近傍で「良い」性質をみたせば, ニュートン法
xn+1=xn−f(xn) f(xn) によってξの近似列{xn}を構成することができる.
● 講義資料
★ 関数へのポインタ 例えば√
2を求めるための区間縮小法のプログラムを考えてみる.
#include <stdio.h>
#define EPSILON 1.0E-10
double interval(double l_side, double r_side, double epsilon) ; int main(int argc, char **argv)
{
printf("%.16f\n", interval(0.0, 2.0, EPSILON)) ; return 0 ;
} /*
*区間縮小を再帰で行う.
* l_side < r_sideの差が epsilonより小さくなったら終了
*そうでなければ,区間縮小を行う
*/
double interval(double l_side, double r_side, double epsilon) {
double c ;
if (r_side - l_side <= epsilon) return (l_side+r_side)/2.0 ; c = (l_side+r_side)/2.0 ; if ((c*c - 2.0) < 0) l_side = c ; else r_side = c ;
return interval(l_side, r_side, epsilon) ; }
ex04.tex,v 1.3 2003-11-02 15:59:52+09 naito Exp
２００３年度・数理解析・計算機数学２・第４回 3
（このプログラムは「再帰」を用いているが,単に「趣味」の問題であるので,気にする必要はない.）この プログラムで,求めるべき真の値αが√2であることを使っているのは,
c*c - 2.0 < 0
の部分のみである.したがって,この部分を独立させてしまえば,
#include <stdio.h>
#define EPSILON 1.0E-10
double interval(double l_side, double r_side, double epsilon) ; double error_sqrt2(double) ;
int main(int argc, char **argv) {
printf("%.16f\n", interval(0.0, 2.0, EPSILON)) ; return 0 ;
} /*
*区間縮小を再帰で行う.
* l_side < r_sideの差がepsilonより小さくなったら終了
*そうでなければ,区間縮小を行う
*/
double interval(double l_side, double r_side, double epsilon) {
double c ;
if (r_side - l_side <= epsilon) return (l_side+r_side)/2.0 ; c = (l_side+r_side)/2.0 ;
if (error_sqrt2(c) < 0) l_side = c ; else r_side = c ;
return interval(l_side, r_side, epsilon) ; }
/* Sqrt[2]がどちらに入るかを計算する*/
double error_sqrt2(double x) {
return x*x - 2.0 ; }
と書き直すことができる.
ここで,interval中で使っている関数error_sqrt2を他の関数に置き換えれば,intervalは他の値を 求めるためにそのまま利用できることがわかる.そこで,関数intervalに引数として「関数」を渡して, 他の値を計算するために利用できないかを考えてみよう.
プログラム中で用いられている関数は,メモリ内のどこかに関数のコードがあり,「Ｃ言語では,関数名
（関数識別子名）はその関数のコードの先頭アドレスを表している」という事実を用いよう.すなわち,「関 数へのポインタ」という変数を定義すれば,それをintervalへの引数として渡すことができる.
「関数へのポインタ」は次のように定義する.「戻り値double,引数はただ一つであり,その型がdouble である関数」へのポインタ変数は
double (*f)(double)
と定義する.これを用いた（上のプログラムの）改良版は以下のようなものである.
ex04.tex,v 1.3 2003-11-02 15:59:52+09 naito Exp
4 ２００３年度・数理解析・計算機数学２・第４回
#include <stdio.h>
#define EPSILON 1.0E-10
double error_cubic_root2(double) ; double error_sqrt2(double) ;
double interval(double, double, double, double (*)(double)) ; int main(int argc, char **argv)
{
printf("%.16f\n", interval(0.0, 2.0, EPSILON, error_sqrt2)) ; printf("%.16f\n", interval(0.0, 2.0, EPSILON, error_cubic_root2)) ; return 0 ;
}
double interval(double l_side, double r_side, double epsilon, double (*f)(double)) {
double c ;
if (r_side - l_side <= epsilon) return (l_side+r_side)/2.0 ; c = (l_side+r_side)/2.0 ; if (f(c) < 0) l_side = c ; else r_side = c ;
return interval(l_side, r_side, epsilon, f) ; }
double error_sqrt2(double x) {
return x*x - 2.0 ; }
double error_cubic_root2(double x) {
return x*x*x - 2.0 ; }
（このプログラムでは,√
2だけでなく21/3も計算している.）
★ ホーナ法
ニュートン法でf(x)が多項式の時には,f(xn),f(xn)を求める必要がある. 多項式f(x) =anxn+ an−1xn−1+· · ·+a0のx=ξでの値f(ξ),f(x)のx=ξでの値f(ξ)を求めるためには,次に述べる
「ホーナ法」がよく用いられる.
多項式f(x) =n
k=0akxkの係数{ak}nk=0から,数列{b(0)k }nk=0,{b(1)k }n−1k=0を次のように与える.
1. b(0)n :=an,
2. b(0)k :=ak−b(0)k+1ξ, k=n−1, . . . ,0.
3. b(1)n−1:=an/n,
4. b(1)k :=ak+1/(k+ 1)−b(1)k+1ξ, k=n−2, . . . ,1.
この時,b(0)0 =f(ξ),b(1)0 =f(ξ)が成り立つ.
これは,f(x)を
f(x) =anxn+· · ·+a0,
f(x) = (x−ξ)(b(0)n xn−1+· · ·+b(0)1 ) +b(0)0 と２通りに書いて,その係数比較をしたものである.
ex04.tex,v 1.3 2003-11-02 15:59:52+09 naito Exp
２００３年度・数理解析・計算機数学２・第４回 5
● 実習内容
•（中）ニュートン法を用いて,√ 2,√
3, 21/3の値を,「相対誤差」10−10で求めるプログラムを,関 数へのポインタを利用して書きなさい.
•（中）区間縮小法を用いて,√n,n1/3の値を,「絶対誤差」10−10で求めるプログラムを,関数への ポインタを利用して書きなさい.（nはプログラム中で与えてよいが,nを変えても本質的にプログ ラムを変更しないようにかきなさい.）
•（易）ニュートン法での(x−1)2(x+ 1) = 0, (x−1)3(x+ 1) = 0の解x= 1への近似の様子をプ ロットしなさい.（特に,f(x),f(x)の計算にホーナ法を用いなさい.）
•電子メールで「今日の講義の感想や意見」を送ってください.
★ 前回の実習の訂正 前回の実習のプログラム
#include <stdio.h>
#define FMT "%02x%02x%02x%02x%02x%02x%02x%02x\n"
int main(int argc, char **argv) {
int i ;
double x=0.0 ;
unsigned char *p ; p = (unsigned char *)&x ;
printf(FMT, *p, *(p+1), *(p+2), *(p+3), *(p+4), *(p+5), *(p+6), *(p+7)) ; for(i=0;i<10;i++) {
x += 0.1 ;
printf(FMT, *p, *(p+1), *(p+2), *(p+3), *(p+4), *(p+5), *(p+6), *(p+7)) ; }
return 0 ; }
は, IntelのCPU上では意図したような表示はしない.
#include <stdio.h>
#define FMT "%02x%02x%02x%02x%02x%02x%02x%02x\n"
int main(int argc, char **argv) {
int i ;
double x=0.0 ;
unsigned char *p ; p = (unsigned char *)&x ;
printf(FMT, *(p+7), *(p+6), *(p+5), *(p+4), *(p+3), *(p+2), *(p+1), *p) ; for(i=0;i<10;i++) {
x += 0.1 ;
printf(FMT, *(p+7), *(p+6), *(p+5), *(p+4), *(p+3), *(p+2), *(p+1), *p) ; }
return 0 ; }
ex04.tex,v 1.3 2003-11-02 15:59:52+09 naito Exp
6 ２００３年度・数理解析・計算機数学２・第４回
が正しいプログラムである.
このプログラムは,６４ビット倍精度浮動小数点数の内部表現をunsigned char型のポインタを用いて, １バイト（８ビット）づつ読み出したものであるが, Intel社（及び, DEC Alpha）のCPUは下位互換性 のため,２バイト以上の数値は「逆順」に格納されている.すなわち,「最上位バイトに数値の下位データ を格納」する.この格納方法を「リトルエンディアン」(little endian)と呼び, IBM, Motorola, Sunなど のCPUは,「最上位バイトに数値の上位データを格納」し,これを「ビッグエンディアン」(big endian) と呼ぶ.1
このように,データの格納ならび順（「バイトオーダリング」(byte ordering)）は, CPUによって異 なるため,数値データのような「バイナリデータ」を交換するためには,バイトオーダリングをチェックす る必要がある.なお,ネットワーク上でデータを交換する場合には「ネットワークバイトオーダリング」と 呼ばれるバイトオーダリングを用いてデータ交換を行うことが規定されている.
なお,Ｃ言語では「ＡＮＳＩ規格内でバイトオーダリングを正確に調べるプログラムを書きなさい」とい う問題は,有名な難問としてよく知られている.
1より正確には,「MSBが下位アドレスにあるもの」がBig Endianであり,「LSBが下位アドレスにあるもの」がLittle Endian である. MIPS社のCPUは, Big/Little Endianを切り替えることができるらしい.
また, “Endian”の語源は,「ガリバー旅行記」に由来する.小人の国では「ゆで玉子を細い方から食べるか,太い方から食べるか」
という議論から戦争が起ったとされている.このとき,「ゆで玉子を細い方から食べる人」を“little end-ian”と呼び,「ゆで玉子を 太い方から食べる人」を“big end-ian”と呼んでいる.
ex04.tex,v 1.3 2003-11-02 15:59:52+09 naito Exp
