２００８年度・数理解析・計算機数学２・第１３回
1
● 前回の講義のまとめ
【２階線型常微分方程式の境界値問題の数値解法】
1.
与えられたp: [0, 1] −→ R, w : [0, 1] −→ R
に対して, x: [0, 1] −→ R
に関する微分方程式x (t) + p(t)x(t) = q(t),
x(0) = 0, x(1) = 0 (1)
を考えよう
.
このように, x
が定義された区間の端点でのx
の条件の下に微分方程式を考える問 題を境界値問題と呼び,
特に,
端点でのx
の値を境界条件として与えたとき, Dirichlet
境界値 問題と呼ぶ.
2.
常微分方程式の境界値問題を解くために,
区間[0, 1]
をN
等分し, t j = jh, h = 1/N
として, {x(t j )}
の近似値を求めることを考えよう.
3.
２階微分はx (t n ) ∼ x (t n+1 ) − x (t n )
h ,
x (t n+1 ) ∼ x(t n+1 ) − x(t n )
h ,
x (t n ) ∼ x(t n ) − x(t n−1 )
h ,
と考えれば
,
x (t n ) ∼ x(t n+1 ) − 2x(t n ) + x(t n−1 ) h 2
によって差分化することができる
. 4.
したがって, (1)
の数値解{x j }
をx j+1 − 2x j + x j−1
h 2 + p(t j )x j = q(t j ), j = 1, . . . , N − 1
によって求めれば良い
.
いま, Dirichlet
境界条件を考慮すれば, x 0 = x N = 0
と置くのが適当で ある.
5.
よって, (1)
の数値解を求めるためには,
連立一次方程式1
h 2 (x 2 − 2x 1 ) + p(t 1 )x 1 = q(t 1 ), 1
h 2 (x j+1 − 2x j + x j−1 ) + p(t j )x j = q(t j ), j = 2, . . . , N − 2, 1
h 2 (−2x N −2 + x N−1 ) + p(t N −1 )x N −1 = q(t N −1 )
を解けばよいことがわかる.
これを行列の形で表現すれば,
A =
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
, P =
p(t 1 )
p(t 2 ) . ..
p(t N −2 )
p(t N −1 )
,
ex13.tex,v 1.1 2008-07-03 17:36:40+09 naito Exp
2
２００８年度・数理解析・計算機数学２・第１３回X =
x 1
.. . x N −1
, Q =
q(t 1 )
.. . q(t N −1 )
( 1
h 2 A + P )X = Q
を解けばよいことがわかる【Gauss-Jordan の消去法】
1. Gauss-Jordan
の消去法とは,
いわゆる「掃き出し法」であり,
連立一次方程式Ax = b
を解く ために,
基本変形行列の積となる行列を利用してMAN (N −1 )x = Mb, MAN =
1 0 · · · 0 0 · · · 0 0 1 · · · 0 0 · · · 0
. . . .. . . . . 0 0 · · · 1 0 · · · 0 0 0 · · · 0 0 · · · 0 . . . . . .. ...
0 0 · · · 0 0 · · · 0
とするアルゴリズムである
.
2. A
が正方行列で正則な場合, Gauss-Jordan
の消去法のアルゴリズムは以下のようなものとなる.
以下でR j は,
その時点での各行列の第j
行目をあらわす.
また, n
は行列のサイズをあらわす.
(a) k = 1, M = I
とする.
(b) (
部分枢軸選択) {a k } n =k たちの中から,
その絶対値が最も大きいものを a k0 としたとき, A, b, M
に対して, R k とR 0 を交換する.
, A, b, M
に対して, R k とR 0 を交換する.
.
(c) (A
の対角成分a kk を1
にする.) A, b, M
に対して, R k ← (a kk ) −1 R k とする.
.
(d) (A
の対角成分以外を0
にする) A, b, M
に対して, = 1
からm
まで, = k
を除いて, R ← R − a k L k を行う.
(e) k += 1
とし,k = n
ならば終了. そうでなければ, Step bに戻る.これが終了した時点で
, b
は解x
を与え, M
はもとのA
の逆行列を与える.
3.
浮動小数点演算による誤差を最小限に止めるために「枢軸選択」は必ず必要な手順である.
4.
これらのGauss-Jordan
の消去法のアルゴリズムは,
行に関する演算については左からA
に基本変形行列を掛けて変形していることに他ならない
.
5. Gauss-Jordan
の消去法で,
正則なn × n
行列を単位行列まで変形するためには, n 2 (n − 1)
回 の乗算とn(n − 1)
回の加算が必要となる.
特にO(n 3 )
回の乗算が必要である.
● 実習内容
• Gauss
の消去法を用いて,
以下の連立方程式を解くプログラムを書きなさい.
3x+12y+9z=3, 2x+ 5y+4z=4,
−x+ 3y+2z=5.
2x+2y+ z=3, x+ y+2z=4, y+ z=5.
ex13.tex,v 1.1 2008-07-03 17:36:40+09 naito Exp
２００８年度・数理解析・計算機数学２・第１３回
3
•
次の行列をLU
分解するプログラムを書きなさい.
また,
その結果を使って上の連立一次方程式を解 くプログラムを書きなさい.
3 12 9
2 5 4
−1 3 2
,
2 2 1 1 1 2 0 0 1
,
• Gauss
の消去法を用いて,
次の常微分方程式の境界値問題の数値解を構成しなさい. (h = 0.01
で解きなさい
)
1. x (t) − sin(πt)x(t) = cos 3 (πt),
x(0) = x(1) = 0.
-0.015 -0.01 -0.005 0 0.005 0.01 0.015
0 0.2 0.4 0.6 0.8 1
x’’ = Sin[Pi*t]*x + Cos^3[Pi*t], H = 0.01
2. x (t) − sin(πt)x(t) = cos 3 (πt),
x (0) = x (1) = 0.
-0.14 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0 0.02
0 0.2 0.4 0.6 0.8 1
x’’ = Sin[Pi*t]*x + Cos^3[Pi*t], H = 0.01
★ 実習内容の解答
