• 検索結果がありません。

● 講義資料

N/A
N/A
Protected

Academic year: 2021

シェア "● 講義資料"

Copied!
5
0
0

読み込み中.... (全文を見る)

全文

(1)

● 講義資料

▼ 講義予定

連立一次方程式の数値解法

掃出し法

– Gaussの消去法

– LU分解

● 前回の講義のまとめ

★ 2階線形常微分方程式の境界値問題

• u: [0,1]−→Rに対する微分方程式の境界値問題

u′′(x) +p(x)u(x) =f(x), x∈(0,1),

u(0) =u(1) = 0 (1)

を離散変数法を用いて数値的に解くことを考える.

定義域の区間[0,1] N 等分して, xk =kh, (h= 1/N)での厳密解の値u(xk)の近似値を uk と書く. この時,2階微分u′′(x)

u′′(x)∼ u(x+h)−2u(x) +u(x−h) h2

と差分化する. これは,

u(x)∼u(x+h)−u(x) h

∼u(x)−u(x−h) h

u′′(x)∼u(x+h)−u(x) h

と考えて導出したものである.

したがって,方程式を差分化すると,

u(xk+1)−2u(xk) +u(xk1)

h2 +p(xk)u(xk) =f(xk) uk+1−2uk+uk1

h2 +pkuk =fk

となる. ここで,この式はk= 1, . . . , N−1 に対して成り立つ式である.

ここで,境界条件u(0) =u(1) = 1は,

u0=uN = 0

(2)

と書き直せるので,解くべき式は u2−2u1

h2 +p1u1=f1, uk+1−2uk+uk1

h2 +pkuk=fk, k= 2, . . . , N−2,

−2uN1+uN2

h2 +pN2uN2=fN2

(2)

となり,未知ベクトルU = (uk)に対する連立一次方程式となる.

いま,

A=

−2 1

1 −2 1

. .. ... ...

1 −2 1

1 −2

, P =

 p1 0

0 p2 0

. .. ... . ..

0 pN2 0 0 pN1

 ,

と書くと, (2)

1 h2A+P

U =F (3)

と書き直すことができる. したがって,2階線形常微分方程式の境界値問題(1) を離散変数 法で解くことは,連立一次方程式(3)を解くことに帰着できた.

以下では,特にp(x) =cの場合を考える. すると,連立一次方程式は 1

h2A+cI

U =F (4)

と書くことができ, 連立一次方程式(4)の可解性は, h2A+cI のすべての固有値が0であ るか否かで判定可能である. いま,Aの固有値(と固有ベクトル)は,

λk=−4 sin2(kπ

2N), Vk=

sin(πkN) sin(2πkN )

... sin((NN2)πk) sin((NN1)πk)

=

sin(jkπN )

である. したがって,h2A+cI の固有値は,N h= 1を使えば, Λk,h=c−4 sin2(kπ

2N) =c−4 sin2(kπh 2 )

と書ける. よって,与えられたh >0に対して, Λk,h6= 0となるとき,任意のF に対して, 立一次方程式(4)の解U を求めることができる.

いま,g(h) =h1sin(h)は,単調減少(hを小さくするほど大きくなる)かつ, limh0g(h) = 1 を使うと,h∈(0, π)において

4 sin2(kπh/2)

h2 >(kπ)2

(3)

が成り立つ. したがって,h∈(0, π)ならば, Λk,h=c−4 sin2(kπh

2 )< c−(kπ)2< c−π2 が成り立つ.

離散変数法で数値解を構成するときには,h >0 は十分小さくとるので, c < π2

が成り立てば, Λk,h<0となり,任意の(十分小さな)h >0に対して(4)は可解となる. に,c <0ならば(4)は可解となる.

★ 連立一次方程式の数値解法

以下,A∈Mn(R),b,x∈Rn として,連立一次方程式 Ax=b

を数値的に解くことを考える. (A が正則でない場合には, 「正則でない」ことが判定でき ればよいことにする.)

連立一次方程式を数値的に解く方法としては,大きく分けて,消去法と反復法がある. はじめ に,消去法の中でももっとも基本的な「掃出し法」(「Gauss-Jordanの消去法」)を考える.

掃出し法とは, 「行基本変形」を繰り返して, 行列 Aを単位行列に変形する方法である. こで,「行基本変形」は,以下の3種類の変形である.

1. i行目をλ倍する.

これは,方程式の両辺に左からDi(λ)をかけることに他ならない.

2. i行目とj 行目を入れ替える.

これは,方程式の両辺に左からTij をかけることに他ならない.

3. i行目にj 行目のλ倍を加える. (ただし,j < iを仮定する)

これは,方程式の両辺に左からEij(λ)をかけることに他ならない.

ただし,行列Di(λ),Tij,Eij(λ)は,

Di(λ) =

i

<

1 0 . . . 0

0 . .. ...

... . .. ...

0 . . . λ . . . 0

... . .. ...

... . .. 0

0 . . . 0 1

< i

(4)

Tij=

i

<

j

<

1 . . . 0 ... . .. . . ... 0 · · · 0 · · · 1 · · · 0

... . .. ...

0 · · · 1 · · · 0 · · · 0 ... . . .. ...

0 . . . 1

< i

< j

Eij(λ) =

j

<

i

<

1 . . . 0 ... . .. . . ... 0 · · · 1 · · · 0 · · · 0

... . .. ...

0 · · · λ · · · 1 · · · 0 ... . . .. ...

0 . . . 1

< j

< i

である.

掃出し法のアルゴリズムは以下の通りである. ただし,「枢軸選択」を行っていないので, 則であっても,アルゴリズムが最後まで到達できないときがある.

– i= 1, . . . , nに対して,以下を繰り返す.

1. 対角成分aii 1となるように,i 行目を定数倍する.

すなわち,「Di(aii1)を左からかける」または,「aij :=aij/aii (j = 1, . . . , n) する」

2. j= 1, . . . , i−1, i+ 1, . . . , nに対して,非対角成分aji 0となるように,j 列目か i列目の定数倍を引く.

すなわち,「Eji(−aji)を左からかける」または,「ajk:=ajk−ajiaik (k= 1, . . . , n) とする」

もし, 1のステップでaii= 0 であるときには,プログラムを停止する.

(5)

● 実習内容

(以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.)

1. (★★★)連立一次方程式





2x+ 2y+ z= 3 x+ y+ 2z= 4 y+ z= 5

Gauss-Jordanの消去法(掃出し法)で解くプログラムを書きなさい.

2. (★★)行列

A=

3 12 9

2 5 4

−1 4 2

の逆行列を求めなさい.

3. (★★★)連立一次方程式





3x+ 12y+ 9z= 3 2x+ 5 y+ 4z= 4

−x+ 3 y+ 2z= 5

Gaussの消去法で解くプログラムを書きなさい.

4. (★★★)連立一次方程式





2x+ 2y+ z= 3 x+ y+ 2z= 4 y+ z= 5

Gaussの消去法で解くプログラムを書きなさい.

5. (★★★)行列

A=

3 12 9

2 5 4

−1 4 2

LU分解を(種々の方法で)求め,それを用いて, b=t(3,4,5),b=t(3,4,−5)に対して, 連立一次方程式Ax=bを解くプログラムを書きなさい.

参照

関連したドキュメント

ニュートン法のアルゴリズム 入力: 関数 とその勾配ベクトル ヘッセ行列 初期点 0 ステップ0: とする ステップ1:

施策前後でツイート内容の傾向に変化があったのかを調べていきましょう。

ここでは, Gauss の消去法で枢軸 選択が必要ない(対角成分が

ここでは, Gauss の消去法で枢軸 選択が必要ない(対角成分が 0 にはならない)場合のみを考える...

ここでは, Gauss の消去法で枢軸 選択が必要ない(対角成分が

この問題では, スツルム列を手で(または, Mathematica などで)計算しておき, 係数を浮動

すなわち, 台形公式と中点則は, 1次以下の多項式に対して正しい積分値を与え,

ここでは, Gauss の消去法で枢軸 選択が必要ない(対角成分が