2005/12/06	数値解析学と演習 福井 先生

先週まで常微分方程式
今週から二回偏微分方程式

# 講議は、来年は一度しかないので、その時にまとめを行う
# 年内は、これと、トピックを一回

==

# 前回言い忘れたので、その内容から ( これは、常微分方程式だけでなく偏微分方程式にも関連する .. )

[陽解法と陰解法]

微分方程式の解法 ( 次の二つある )
	陽解法
		# 例 : ルンゲ＝クッタ / オイラー法
			y(x+h)=y(x)+h f(x,y(x)))
			この式の右辺は、x の点の値のみを利用している
				y(x)
				h
				f(x,y(x))
					どれも、x の時の値
			# 解っている値だけを利用して、次の計算を行う
			## 値の計算では、未知の値が存在しない

	陰解法
		# 例 ( 説明のために作った )
			y(x+h)=y(x+h)+h f(x,y(x)))
			       ^^^^^	<= 右辺に解っていない値を利用する公式
			# 「理論上」は、値の計算の右辺に、解っていない値を利用することができる。
			## 計算は、「直接は」できない ( 未知の値が入っているので.. )
				=> 方程式を解く必要がある
				# 一般には、線型化し、連立方程式を解くという手法を取る

		# Text には、ルンゲ=クッタの陰解法の例が乗っている (text p.528)

		text p.552 ギア法
			各々のステップでの計算で、未知の値があるので、連立方程式を解く必要がある
				# これは、今日、学ぶ予定の偏微分方程式でも同様
			# ギア法は、これだけの本があり、そこにプログラムも乗っている

	[陽解法と陰解法の性質の違い]

		陽解法	簡単			Δt (独立変数) を大きく取れない
						# 大きくすると不安定になる
						## 偏微分方程式で例を示す
						誤差が積み重なる

		陰解法	方程式を解く必要がある	Δt (独立変数) を大きくできる
						安定性がある


		例 : 硬い系 ( スティッフな系 )
			時定数が大きく異る系
				低周波と高周波の和 ( 三桁位異る場合 )

			周期関数を近似するには、ある程度の細かさで刻む必要がある
				# 最低限 Δt < \frac{ω}{10}
				## 周期が一つならば、これで Okey
				### 二つの周波数があると、小さい方に合せる必要がある
				###	=> Δt を、あまり大きくできない
				###	しかし、Δt を小くすると、「有限桁」の問題が..

	# 陽解法で、Δt を小くした結果の計算量の増大と、陰解法の方程式を解く計算量の増大は、同程度かもしれない

	## 後で、偏微分方程式の例では、陽解法が不安定になる例をやる

==

[偏微分方程式]

	常微分方程式と偏微分方程式の違い
		=> 独立変数が二個以上ある

		# 物理現象では、最大 4 つ ( 空間の三次元 x,y,z と時間 t )
		## 数学的には、何次元でも..

	例
		ラプラス方程式
		\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}

	[2 階線型偏微分方程式]

# 以下、
# \PD{u}{x} = \frac{\partial u}{\partial x}
# \PDD{u}{x} = \frac{\partial^2 u}{\partial x^2}
# \PDC{u}{x}{y} = \frac{\partial^2 u^2}{\partial x\partial y}

	a\PDD{u}{x} + b\PDC{u}{x}{y} + c\PDD{u}{y} + d\PD{u}{x} + e\PD{u}{y} + f = 0

	頭の三つの項目によって、分類される
		D = b^2 - 4ac
	とした時
		D<0 楕円型
			\PDD{u}{x}+\PDD{u}{y} = 0	ラプラス方程式
			\PDD{u}{x}+\PDD{u}{y} = f	ポワソン方程式

			# 太鼓の皮の振動

			## 3 次元であれば、
			##	\PDD{u}{x}+\PDD{u}{y}+\PDD{u}{z} = 0

		D=0 放物線型
			\PD{u}{t} = \PDD{u}{x}		熱方程式

		D>0 双曲型
			\PDD{u}{t} = \PDD{u}{x}		波動方程式

			# 海の波の方程式

[偏微分方程式を解くための条件]
	初期値	t (時間) に対する条件
		[例]
			箱の一方向が空いている
				空いている辺は自由端になる
					最初の t の時にどうのような状態かを示す
						=> 初期状態

	境界値	x,y,z (空間) に対する条件
		[例]
			蓋のない箱にゴムを被せる
				箱に接している部分は動かない => 境界条件

	# 条件が少いと、自由すぎて、答がでない
	#	境界値が現象を決めている
	#		方程式の系が違っていても、境界条件が同じだと答が似る
	#		方程式の系が同じでも、境界条件が違うと答ががらりと変る

	## 一点しか固定されていない紐は自由に動ける
	##	二点を固定すると答の取り得る範囲は限られている !!

	[値の与え方]
		1) (関数)値 u を与える
			第1種境界条件 ( ディレクレ : Dirichlet )
			# 境界の値を 0 にする

		2) 法線方向の微係数 \PD{u}{n}
			第2種境界条件 ( ノイマン : Neumann )
			# 境界で、熱が逃げる速度を示す ( 温度差のα倍 )
			## 熱方程式などはこちらが自然

		3) u と \PD{u}{n} の 1 次式
			第3種境界条件
		

	ディレクレ条件の場合は、誤差が出にくい
		境界では誤差がまったくないので、計算するにつれて誤差が補正される
		# 途中で、多少の誤差があっても、境界に影響されて誤差が失われる

	ノイマン条件の場合は、自由に移動できる端が出る
		誤差が累積する傾向がある

==

	偏微分方程式を解くには、境界条件が必須
		海底の形状などは、情報が少い ( 地上ならば、国土地理院が頑張っているが.. )


	# 特別な例
	#	境界条件のない偏微分方程式の例
	#		社会現象を記述する場合
	#			境界条件が縮退するような仕組で計算する
	#				# 物理現象とは異る挙動になっている

	波のシミュレーション
		東京湾や大阪湾の模型があり、実際に実験を行う

	# 偏微分方程式の解は、計算結果が本当に正いかどうかが解らない
	#	測定結果と比較するしかない
	#	[例] 淡路島の海流の計算結果
	#		=> 昔の海軍の測定結果位しかないので、比較するのが大変

==

[解法] ( 解析的に解けるのが一番なのだが、そうなることは稀なので.. )
	解析的に解く ( 解ける場合も、形状が複雑だと駄目 )
		# 偏微分方程式は、方程式の型だけでなく、問題の形状に大きく影響される
	数値的に解く ( 離散化を行う )
		(有限)差分法 (FDM)
			微分を差分に置き換えて解く

			[例] ポワソン方程式 \PDD{u}{x}+ \PDD{u}{y}=0 を考える
				# これを、二階の中心差分で、差分化すると..

			\frac{u(x+Δx,y)-2u(x,y)+u(x-Δx,y)}{Δx^2} +
			\frac{u(x,y+Δy)-2u(x,y)+u(x,y-Δy)}{Δy^2}

			となるので、u(x,y) を u(x+Δx,y) 等の値から計算する

			係数行列の内容は、3重対角+2斜線行列になる ( この場合は二次元なので..)
				# 境界値は、境界条件によって定まる
				対角と対角の上下は零にならない
					=> スパースな行列になる

			# 有限要素法で、差分の方法と分割方法が同じなら、同じ係数行列になる。

			## 三次元だと、3 重対角+4斜線行列になっている


		有限要素法 (FEM)
			問題(領域)を、小さな要素(矩形や三角形)に分解して解く
			その要素の中での関数値を近似する ( 例 : テント関数 )
				=> 要素間の関係が連立方程式(行列)の形で定式化される。
					方程式を解くと答が出る
			# ポイントは、「要素分割」の方法

			[手順]
				1) 要素の分割
				2) 各々要素に答(の形)を仮定する
					通常は多項式 ( 下手をすると線型 )
					# sin/cos でも良いが計算が大変なので..
					係数を未知数として扱う
				3) 各々の要素間の関係式(係数行列)を作る
					# 個々の要素の答がきまるので、係数だけを全体の条件に合せて、求める
					係数行列が境界条件を満すようにする
				4) 連立一次方程式を解く
		境界要素法 ( 電気の分野位でしか使わない )
			# 電気の導体の中の電位を計算するために考えられた
			## => 代用電荷法 ( これの一般化したものが境界要素法 )

	解法による行列の形の違い
		(有限)差分法 (FDM)
		有限要素法 (FEM)
			=> 疎(スパース)行列になる
		境界要素法
			=> 密(デンス)行列になる


# 解法としては、それほど難しい話はないが、実際にプログラムを作成すると
# なると、話は別で、仔魔界所が色々と大変

==

今、特定な係数行列を与えると、どの解法が望ましいかを答えてくれる web site を作っている所

==

[演習] (福井先生が解く)
	# プログラムを作ると一般に大変だが、この例は、Excel でできる
		# Δt を大きくすると、不安定現象がおきる ( 負の数が出る )
		## ぎざぎさになる

	# 一次元なら、Excel で十分、二次元でも、なんとか..
	## 式が複雑過ぎるとだめだけど..
	### オイラー法も、Excel で計算できる ( 一列だけなので、あまり面白くない )

	# 「Excel で色々な問題を解く方法」という本があるので興味があれば..

	問題は、
		\PD{u}{t}=a^2\PDD{u}{x}
	を考える。

	# 差分化すると..

		\PD{u}{t} = \frac{u(t+Δt,x)-u(t,x)}{Δt}
		\PDD{u}{x} = \frac{u(t+Δt,x)-2u(t,x)+u(t-Δt,x)}{Δt^2}

	これを、上の式に代入すると..

\frac{u(t+Δt,x)-u(t,x)}{Δt} = \frac{u(t,x+Δx)-2u(t,x)+u(t,x-Δx)}{Δx^2}

	これを変形して、u(t+Δt,x) を求める式に変形すると

u(t+Δt,x) = u(t,x) + \frac{Δt}{Δt^2}(u(t+Δt,x)-2u(t,x)+u(t-Δt,x))

	この式は u(0,x) が初期値条件として与えられていれば、左辺を求めるための右辺の値は、前以て計算できるので、これは、陽公式。よって陽解法によって計算できる

	初期値は、sin で与えてある。境界条件は 0 ( ディレクレ条件 )

==

# 本当に、偏微分方程式を解くだと、色々なことを知らないとだめ ( 講議は、入門程度.. )

偏微分方程式を解く方法 (三つ)
	差分法/有限要素法/境界要素法 ( 普通は、最初の二つ )
		連立方程式が陽解法なら簡単だが、陰解法の場合は、連立方程式を解く必要がある

# 差分法と有限要素法の違いは、微分(強微分/強形式)と、弱微分/弱形式の違い

## 有限要素法は、元々は建築関係から出た計算法だが、後に数学的に定式化された

一次元の問題を差分化すると、行列は、三重対角になる
ニ次元の問題を差分化すると、行列は、三重対角+二斜線(位置は、分割の方法で決る)
	できるだけ、要素が真中にあつまっている方が計算が楽
		# 長方形の分割の場合は、数の少い方向から番号付けすれば、真中に集る

	陽解法と陰解法
		陽解法の場合は、安定条件を満すような小さな h しか取れない