講義ノート

第 1 回: 数値計算の基礎 (10.06.2010)

数値解析 (Numerical Analysis)

1. 数値計算
2. 関数近似
3. 数値積分
4. 線形方程式
5. 行列の固有値問題
6. 非線形方程式

参考図書:

1. 数値計算の基礎と応用
2. 数値計算の基礎
3. だれにでもわかる数値解析入門
4. Maxima

計算コスト:

1. 計算時間
2. 使用メモリ量
3. 近似精度 (連続量の場合)

1. 精度と誤差
2. 有効桁数
3. 計算機上の浮動小数点数
4. 正規化表現の限界
5. 丸めと丸め誤差
6. 単精度と倍精度の比較
7. 関数による誤差伝播 (1 変数関数) テーラー展開
8. 関数による誤差伝播 (多変数関数) テーラー展開
9. 四則演算の誤差伝播 (加減算, 乗算, 除算)

第 2 回: 関数計算(展開と補間) (10.13.2010)

関数を近似しない

たたみこみ

ラグランジュ

ニュートン補間

ネスティング法

問 4 (プログラミングをすること, 発表)

第 3 回: 数値積分

第 4 回: 線形変換の誤差解析

ノルムが小さいこととすべての成分が小さいことは等価

1-ノルム (証明をみること)

2-ノルム

OO-ノルム

絶対誤差と相対誤差はスカラーである

絶対誤差が小さいことは, X(ベクトル) が X(ベクトル) に成分ごとに近いことと等価

ノルムの公理(距離の公理)

公理: 根幹的なもの

定理: 証明されべきもの

定義: ルールを定義する

系: 証明されべきもの

補題: 証明されべきもの (定理証明ようのいくつかの小さい問題)

ノルムの等価性

ベクトル列 (lim に参照すること)

ノルムの増幅率

ノルムがスカラーである

ベクトルノルムと行列のノルムを混同しないようにちゃんと区別しよう

従属ノルムの性質 (1, 2)

逆変数を持つ線形変換 A の最小増幅率

増幅率の最小値

最小増幅率は

A が対称行列 (?)

直行行列およびそれに 0 でないスカラー倍した行列は, 2-ノルム

問題 4

相対誤差を示すこと

有効桁数を示すこと

第 5 回: 線形方程式の直接解法

第 6 回: 線形逆変換の構成 (2010.12.08)

演習: 問 3 (趙飛)

ピボット選択の効果 (誤差を小さくするために)

LU 分解 (A行列だけをみる, b ベクトルを気にせずに)

1. ピボット選択
2. 正則判定
3. 行並べ替える (行列成分を 0 になるために)
4. r を替える
5. 出力 A p (ピボット)

前回より手間かかるが,

単位下三角行列 (対角成分は 1 となる)

Ax = b = Lc = LUx

任意のベクトル x について成立するので, A = LU

行列の下三角行列と上三角行列の積への分解 => LU 分解 (計算が楽になる)

列対角優位

行対角優位

余計なものが入っているために, 注意する

コレスキー分解法

第 7 回: 帯係数行列の線形方程式

第 8 回: 最小二乗法 (2011.01.05)

最小二乗法: 連立方程式の残差の 2 -ノルムを最小にする解を求める方法

幾何学における垂線の問題

実験データの直線近似, 曲線近似.

(Leaset Squares Method)

実験データなどの誤差を含んだ値から, 最もフィットする関数を計算する手法が最小二乗法 (method of least squares) です. 一番簡単なのは, 直線で近似する方法です.

y = ax + b

a = ?

b = ?

数値計算ではこれらの式から a と b の値を計算する. なぜ [最小二乗] なのかというと, 上の式を求めるときに [直線と各データポイントの差の二乗を最小にする] という手法が使われているからです.

GSL (GNU Scientific Library) はその名のとおり, 科学技術計算のライブラリで, さまざまな数値計算法の関数がたくさん集めれている. ANSI C で記述されていて, C や C++ から呼び出せる. 次のような計算法ものがサポートされており, 数値計算の教科書に載っているようなものは大体使える.

複素数 (Complex Numbers) 多項式 (Polynomials) 特殊関数 (Special Functions) ベクトルと行列 (Vectors and Matrices) 順列 (Permutations) 組み合わせ (Combinations) ソーティング (Sorting) BLAS のサポート (BLAS Support) 線形代数 (Linear Algebra) 固有空間 (Eigensystems) 高速フーリエ変換 (Fast Fourier Transforms, FFTs) 数値積分 (Numerical Integration) 乱数の生成 (Random Number Genertion) 準ランダム数列 (Quasi-Random Sequences) ランダム分布 (Random Number Distributions) 統計 (Statistics) ヒストグラム (Histograms) モンテカルロ法 (Monte Carlo Integration) N タプル (N-Tuples) 類似アニーリング (Simulated Annealing) 常微分方程式 (Ordinary Differential Equations, ODE) 補間法 (Interpolation) 数値微分 (Numerical Differentiation) チェビシェフ近似 (Chebyshev Approximations) 級数加速 (Series Acceleration) 離散ハンケル変換 (Discrete Hankel Transforms) 一次元の根の決定 (One dimensional Root-Finding) 一次元の極小化 (One dimensional Minimization) 多次元の根の決定 (Multidimensional Root-Finding) 最小二乗法 (Least-Squares Fitting) 非線形最小二乗法 (Nonlinear Least-Squares Fitting) 物理定数 (Physical Constants) IEEE 浮動小数点計算 (IEEE Floating-Point arithmetic)

線形システム Ax = b

数式処理システム (Maxima) 数値計算システム (MATLAB, Octave, RLab, Scilab) 統計処理 (LISP-STAT, R) 幾何関係 (Cinderella, KSEG, dynagraph, Geomview, SnapPea, Surf, surfer, Teddy, )

変数と関数の定義:

変数名: 変数の内容

kill (値を消去する変数名 1, 値を消去する変数名2, ...);

関数名 (変数名):=関数式;

Cobb-Douglas Production Function (コブ=ダグラス型生産関数)

solve(x^3+2*x^2+3*x+4=0,x);

これをニュートン法を用いて, 数値的に解を近似する方法がある. そのために, 必要な計算パッケージ newton.mac を読み込ませる.

パッケージを読み込むには

load(パッケージ名);

mtx:matrix([2,-1,0,0],[-1,2,-1,0],[0,-1,2,-1]); transpose(mtx);

Maxima プログラミングの初歩(基本的な手順は次の通り):

1. 新しいユーザ関数は, 何を入力して, 結果として何を出力するかを明確にする
2. その結果を得るために, 手作業ではどのような手順を踏むか, 細かくて丁寧に考える.
3. ある処理の結果をそのまま次の手順の入力データとするならば, ネストさせる.
4. 複数の処理の結果を引数とするならば, 各処理結果を変数に保管しておく.
5. 一通り完成したら, いくつかの値を代入してテストする. 特に, [いじわる] な値を選ぶといいでしょう.

fpprec:30; bfloat(sqrt(%pi/2));

log(%e); -> logarithmic function

sum(i, i, 1, 10) sum(式, 変数, 開始値, 終了値);

固有値 eigenvalues

制約条件がない最適化

plot2d([0,x^6-6*x^4+13*x^2-2],[x,-2,2]);

非線形方程式を解くためのアルゴリズム:

1. 2 分法 (Bisection Method)
  長所: 閉空間 [a, b] に解があれば, 必ず, 解に収束する. 間違いなく解を探すので, ロバストな解法と言われている.連続であればどんな形の関数でも解に収束するので信頼性が高いのである. さらに, 解の精度も分かり便利である.
  短所: 収束が遅い. (一次収束)

2. はさみうち法(False Position Methond)
3. 割線法 (Secant Method)
3. ニュートン法 (Newton Method)
  長所: 初期値が適当ならば, 収束が非常に早い. (二次収束)
  短所: 初期値が悪いと, 収束しない. 収束しない場合がるので, 反復回数の上限を決めておく必要がある.

2 月 9 日 (水, 12:00) 出さないと, 22 点をもらえる

最後のレポートに

test1.bat ファイルを作る test2.bat ファイルを作る

batch(“test1.bat”);

remvalue(all)$ remfunction(all)$ load(numericalio.lisp)$

f(x):=sin(x)+cos(x)$ x[0]:-1.0$ x[1]:2.0$

eps:1.0e-4$ k:1; dx(x):=-f(x[k])/(f(x[k]) -f(x[k-1]))$ for i : 1 while dx(x[k])*dx(x[k-1]) > eps do (if k < 10000 tehen (del : ))

batch(“test2.bat”);

1.0e-4

work.mac htop slidy pyrect

newton_methond():=

block(

for i:0 do (

)

);

batch(“test.mac”)

第 9 回: 非線形方程式

第 10 回: 数式処理ソフト利用法