(4) 14 回微分方程式の数値解法 電 301 数値解析第

電 301 数値解析 第 14 回

微分方程式の 数値解法 (4)

電301数値解析(2017)琉球大学工学部電気電子工学科 担当:半塲 1

5

(1)

•

2

v

+ v

+ v

+ v

− 4v

= h

f

を紹介したが, これを行列で表現する方法の 一般はまだ説明していなかった.

•

5

(2)

格子点 5点差分公式

g₇ g₈ g₅ v₂₁ v₂₂ g₆ g₃ v₁₁ v₁₂ g₄

g₁ g₂

g1+v12+v21+g3−4v11=h²f11

g₂+g₄+v₂₂+v₁₁−4v₁₂=h²f₁₂ v₁₁+v₂₂+g₇+g₅−4v₂₁=h²f₂₁ v₁₂+g₆+g₈+v₂₁−4v₂₂=h²f₂₂

5

(3)

格子点 5点差分公式

v₂₁ g3 v11 v12

g₁

g₁+v₁₂+v₂₁+g₃−4v₁₁=h²f₁₁ 第1行目の公式では左の部分を参照する. v₁₁のまわりの要素を足し合わせてから v11の4倍を引いたものがh²f11

と等しいという式. 第2行目以降も同様.

5

(4)

この場合の5点差分公式の行列表現は以下の通り.











4 −1 −1 0

−1 4 0 −1

−1 0 4 −1

0 −1 −1 4



















 v₁₁ v₁₂ v21

v₂₂











=









 g₁+g₃ g₂+g₄ g5+g7

g₆+g₈











−h²









 f₁₁ f₁₂ f21

f₂₂











vijがv₁₁,v₁₂,v₂₁,v₂₂という順で並べられていることが重要. これを一般化する([山本]).

5

(5)

• x軸,y軸にそれぞれNx個,Ny個の格子点が取られているものと する. 前回と同様に,格子点を, (1,1), . . . , (Nx,1), . . . , (1, Ny), . . . , (Nx, Ny)であらわす.

• 境界を, (1,0), . . . , (Nx,0), (1, Ny+ 1), . . . , (Nx, Ny+ 1) (0,1), . . . , (0, Ny), (Nx+ 1,1), . . . , (Nx+ 1, Ny)とし(次ページ図),縦 ベクトルv= (v1, . . . , vNxNy)^Tをその次のページのように定義 する. 続いて,同じ並べ方で縦ベクトルfを構成する. さらに, 行列A,ベクトルgを後述のように構成すると・・・

5

Av = g − h

f

(1,1) (2,1) (Nx,1) (1,2) (2,2)

(1,Ny) (2,Ny) (Nx,Ny)

(Nx,2)

(1,0) (2,0) (Nx,0)

(1,Ny+1) (2,Ny+1) (Nx,Ny+1)

(0,1) (0.2) (0,Ny)

(Nx+1,1) (Nx+1,Ny)

(Nx+1,2)

U

(1,Ny+1) (2,Ny+1) (Nx,Ny+1)

(0,1) (0.2) (0,Ny)

(Nx+1,1) (Nx+1,Ny)

(Nx+1,2)

v

v

v

v

v

y-1)Nx+1

v

yNx

Nx次の正方行列Hを次式左のように定義する. また,IをNx次の 単位行列とする.これらを使って,NxNy次の正方行列Aを次式右の ように定義する.

H=















 4 −1

−1 4 −1 . .. ... ...

. .. ... −1

−1 4















 ,A=















 H −I

−I H −I . .. ... ...

. .. ... −I

−I H

















g= (g01+g10, g20, . . . , gNx−1,0, gNx,0+gNx+1,1;g0,2,0, . . . ,0, gNx+1,2; . . .;g0,Ny−1,0, . . . ,0, gNx+1,Ny−1;

g0,Ny+g1,Ny+1, g2,Ny+1, . . . , gNx−1,Ny+1, gNx,Ny+1+gNx+1,Ny)^T とする.

Lax

(1)

•

(放物型,

•

([河村]).

Lax

(2)

•

∂u

∂t = Lu

(x, t)

L

(x

Lax

(3)

•

u(x

, t

)

v(x

, t

)

(t

− t

= ∆

, x

− x

= ∆

).

∆xと∆tは独立ではなく,ある関数h(ただしlimθ→0h(θ) = 0) が存在し, ∆x=h(∆t)となっているものとする.

• ∂u

∂t = Lu

v(x

, t

) = S(∆

, ∆

)v(x

, t

)

Lax

(4)

• u(x, t + ∆

) − S(h(∆

), ∆

)u(x, t) = o(∆

)

S(∆x,∆t)はLに対応する差分演算子である.

f(θ) =o(θ)という記法は, limθ→0,θ6=0f(θ)/θ= 0を意味する.

• ∃T > 0, ∀n, ∀∆

, 0 < n∆

< T ⇒ kSk < C

S

kSk

Lax

(5)

• Lax

v (x

, t

) = S(∆

, ∆

)v (x

, t

)

∂u

∂t = Lu

→ 0

v(x

, t

)

([河村]).

放物型偏微分方程式の数値解法

(1)

•

1

(x, t)

(x ∈ [0, 1], t ∈ [0, T ]),

2

t = 0

u(x, 0) = a(x)

u(0, t) = 0,

u(1, t) = 0

•

[田端][河村]

放物型偏微分方程式の数値解法

(2)

∂u

∂t = ∂

u

∂

x +f, u(x, 0) = a(x), u(0, t) = u(1, t) = 0.

x=0

t

u(x,0)=a(x)

x=0

x=1 x=1

t=0 t=T

• 2変数の定係数2階線形放物型偏微分方程式は^∂u_∂t =a^∂_∂²2^ux+f であるが, τ = atとすると, ^∂u_∂τ = ^∂u_∂τ^∂τ_∂t = a^∂u_∂τ であるから, a^∂u_∂t =a^∂_∂²2^ux+fとなり,この両辺をaで割ってf¯=f /aとすれ ば, ^∂u_∂t =^∂_∂²2^ux+ ¯fという形になる. この講義ではこの形への変 形はすでに済んでいることを前提としている.

• 熱伝導方程式の場合は,f(x, t)は外部熱源に相当する.

• 境界条件がu(0, t) =g0(t),u(1, t) =g1(t)となっている場合には,

¯

a(x) =a(x)−(1−x)g0(0)−xg1(0), ¯f=f−(1−x)g0^′(t)−xg1^′(t)と して,^∂w_∂t =^∂_∂²2^wx+ ¯f, w(x,0) = ¯a(x), w(0, t) =u(1, t) = 0を解き, u(x, t) =w(x, t) + (1−x)g0(t) +xg1(t)とすれば,もとの問題の 解が得られる. よって, ∂u

∂t =∂²u

∂²x+f, u(x,0) =a(x), u(0, t) =

u(1, t) = 0.という単純化した問題を考えても一般性は失われ

ない.

放物型偏微分方程式の数値解法

(4)

•

=

+ f

•

Euler

•

/(∆

)

≤ 1/2

放物型偏微分方程式の数値解法

(4)

•

Euler

=

∂²v

∂²x

+f

, t

)

(x

, t

)

θ

θ

0

1

v(x

, t

) − v(x

, t

)

∆

= θ v(x

, t

) − 2v (x

, t

) + v(x

, t

)

∆

+ (1 − θ) v(x

, t

) − 2v(x

, t

) + v(x

, t

)

∆

+ θf (x

, t

) + (1 − θ)f(x

, t

)

放物型偏微分方程式の数値解法

(6)

• θ = 0

Euler

θ = 1

Euler

θ = 1/2

Crank-Nikolson

• θ ≥ 1/2

∆

/(∆

)

≤

1/2

θ = 1/2

Crank-Nikolson

([河村]).

双曲型偏微分方程式の数値解法

(1)

• 2

2

= c

•

(

− c

)(

+ c

)u = 0,

(

+ c

)(

− c

)u = 0,

双曲型偏微分方程式の数値解法

(2)

•

− c

= 0,

+ c

= 0

2

•

g

(x − ct)

g

(x + ct)

g

(·)

(i = 1, 2).

双曲型偏微分方程式の数値解法

(3)

•

− c

= 0

+ c

= 0

•

x

v(xj,tn+1)−v(xj,tn)

∆t

+ c

x

= 0.

双曲型偏微分方程式の数値解法

(4)

• µ = c∆

/∆

v(x

, t

) = (1 − µ)v(x

, t

) + µv(x

, t

)

µ

(Courant Num- ber)

•

µ ≤ 1

(ただし, ∆

> 0, ∆

> 0

双曲型偏微分方程式の数値解法

(5)

•

µ > 0

µ = 0

v(x

, t

) = v (x

, t

)

•

µ = 1

(これは許容される),

の公式は

v(x

, t

) = v(x

, t

)

(g

(x − ct))

差分法に関する補足

(1)

•

•

的

Euler

があることが知られている

([河村]).

差分法に関する補足

(2)

•

差分法に関する補足

(3)

•

([河村])

有限要素法

(1)

•

• [Gupta and Meek]

Courant(1942)

Ar- gyris (1954), Turner (1956), Clough (1960), Zienkiewicz and Cheung (1965)

上記で

Courant

有限要素法

(2)

•

•

x = (x, y)

x = (x, y)

有限要素法

(3)

•

R

D

uv dx

(u, v)

u

v

R

D

u

v d x

( u , v )

• u

∇u

有限要素法

(4)

• ∂

u

∂x

+ ∂

u

∂y

= f (x, y ) ((x, y) ∈ U ), u = 0 ((x, y) ∈ ∂U)

a

= a

= 1, a

= a

= 0, b

, b

= 0, c = 0

•

∇v) = (f, v)

有限要素法

(5)

•

•

∂U

有限要素法

(6)

•

∀v, (∇u, ∇v) = (f, v)

u

u, v

{φ

: 1 ≤ i ≤ N }

(u = P

i

p

φ

, v = P

i

q

φ

{c

: 1 ≤ i ≤ N }

Galerkin

有限要素法

(7)

• Galerkin

u, v

P

i,j=1

p

q

(∇φ

, ∇φ

) = P

j=1

q

(f, φ

)

a

= (∇φ

, ∇φ

),

b

= (f, φ

)

A = (a

)

a

b , p , q

p

Aq = b

q

A

有限要素法

(8)

• v(の展開係数)

∀q, p

Aq = b

q

p

Ap−b = 0

有限要素法

(9)

•

φ

•

U

• R

R

有限要素法

(10)

R

(とくに三角

U

U

有限要素法

(11)

•

φ

(要素が

•

3

1

1

2

3

3

x y u

1

x y u

1

x y u

1

要素 要素 要素

このようにすると, 隣接する要素は節点を共有し, 節点で

1

隣接要素で 内挿関数を作り まとめる

有限要素法

(14)

• U

∂U

• ∂U

(∂U

u = 0)

有限要素法

(15)

•

N

≤ i ≤ N

i

•

i

K

Γ

, . . . , Γ

有限要素法

(16)

• 1 ≤ k ≤ K

k

i

1

U \ ∪

Γ

φ

φ

有限要素法

(17)

• 節点として要素の頂点を選択しないこともある. この 場合,隣接する要素Aと要素Bの共有する頂点をcと したとき,要素Aで定義され頂点cで零とならない内 挿関数φAと要素Bで定義され頂点cで零とならない 内挿関数φBが連続に結合できないことがあり得る. こ のような内挿関数の取り方を非適合という. 非適合な 内挿関数が利用されることもある.

有限要素法

(18)

•

Galerkin

a

= (∇φ

, ∇φ

), b

= (f, φ

)

Ap − b = 0

p = (p

, . . . , p

)

v = P

i=1

p

φ

•

(文献を参照).

有限要素法

(19)

• φ

Γ

, . . . , Γ

≤ k ≤ K

∇φ

, ∇φ

)

• (f, φ

)

性がある.

有限要素法

(20)

• 後で必要になるので, 1次元の有限要素法における要素, 節点,内挿関数,φiについて述べる. ただし,内挿関数 は線形関数とし,適合の条件を満たす内挿関数に議論 を限定する.

• 1次元の有限要素法では,領域はRⁿの(有界)閉区間で, 節点はその区間に取られた有限個の点である. 節点が x₁, . . . , xNという順に並んでいるものとすると,要素は 閉区間[xi, xi+1]である(1≤i≤N−1).

有限要素法

(21)

• 内挿関数としては,各要素ごとに, その左端で1, 右端 で0となる関数と,その左端で0,右端で1となる関数 を取る.

• 各節点で, その点で1となる内挿関数を集めて関数φi

を作ることは2次元の場合と同じ.

要素1 要素2 要素3 節点1 節点2 節点3 節点4

0

u 1

関数1-2

関数2-1

関数2-2 関数3-1 関数3-1

節点1 節点2 節点3 節点4

有限要素法

(23)

•

[Quarteroni].

1

a

= 1, b

= 0, c = 0

∂u

∂t , v

+ ∂u

∂x , ∂v

∂x

= (f, v)

有限要素法

(24)

•

N

v

{φ

(x) : 1 ≤ i ≤ N }

(x

1

φ

(x)

v = P

i=1

q

φ

(x)

• u =

N

X

i=1

u

(t)φ

(x).

(係数が時変).

有限要素法

(25)

• m

= (φ

, φ

), a

= (

,

), b

= (f, φ

)

M = (m

)

, A = (a

)

を, これらの要素をまとめた行列とする. ま た,

b = (b

, . . . , b

)

u(t) = (u

(t), . . . , u

(t))

, q = (q

, . . . , q

)

•

u, v

有限要素法

(26)

• ( du

dt )

M q + u

Aq = b

q.

q

M , A

M d u

dt + Au = b

ち, 偏微分方程式が常微分方程式に帰着され たことになる.

有限要素法

(27)

上述の常微分方程式の数値解法は何でもよいので あるが, ここでは

θ

θ

M u(t

) − u(t

)

∆

+ A (θu(t

) + (1 − θ))u(t

))

= (θ b (t

) + (1 − θ) b (t

))

有限要素法

(28)

•

θ

t

t

• θ

θ = 0

Euler

θ = 1

Euler

θ = 1/2

Crank-Nikolson

有限要素法

(29)

•

θ ≥ 1/2

≤ θ < 1/2

Aw = λM w

([Quarteroni]

有限要素法

(30)

•

[Grossmann et al].

a

= 1, b

= 0, c = 0

∂

u

∂t

, v

+ ∂u

∂x , ∂v

∂x

= (f, v)

有限要素法

(31)

•

u =

N

X

i=1

u

(t)φ

(x).

する

(係数が時変).

とベクトルを用いることにすると, 解くべき 方程式は,

M d

u

d

t + Au = b

有限要素法

(32)

•

1

(連立)

• ∂u

∂t +a ∂u

∂x +bu = f

v

, v

+ a

, v

+ b(u, v ) = (f, v)

(上述の偏

有限要素法

(33)

•

([Quarteroni]).

streamline diffusion finite element method,

不連続

Galerkin

といった手法が用いられる

([Quarteroni]).

有限要素法

(34)

•

Galerkin

([Li and Chen]).