偏微分方程式の数値解析
お忙しいところすみません.
現在,FDTD法による数値解析を行っています.
FDTD法は空間及び時間微分に対して中心差分を行っていますがこれを4次のルンゲクッタ法で行うことは出来るのでしょうか?ある参考書に以下のような「連立1階微分方程式」を4次のルンゲ-クッタ法で解くというのがありました.
dy/dt=f(t,y,z)
dz/dt=g(t,y,z)
これを4次のルンゲ-クッタ法で解くと
y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4)
z(n+1)=z(n)+1/6*(j1+2*j2+2*j3+j4)
ここで
k1=Δt*f(t(n),y(n),z(n))
j1=Δt*g(t(n),y(n),z(n))
k2=Δt*f(t(n)+Δt/2,y(n))+k1/2,z(n)+j1/2)
j2=Δt*g(t(n)+Δt/2,y(n))+k1/2,z(n)+j1/2)
k3=Δt*f(t(n)+Δt/2,y(n))+k2/2,z(n)+j2/2)
j3=Δt*g(t(n)+Δt/2,y(n))+k2/2,z(n)+j2/2)
k4=Δt*f(t(n)+Δt,y(n)+k3,z(n)+j3)
j4=Δt*g(t(n)+Δt,y(n)+k3,z(n)+j3)
これをマクスウェルの方程式に対応させ
ε∂E(t,x,y,z)/∂t=f(t,H)=∇×H
μ∂H/∂t=g(t,E)=-∇×E
として出来ないこともない気がするのですが.
なぜこの方法を試みようとした理由は現在,通常のFDTD定式化により数値解析を行っていますが
・プログラムは正しい
・始めの部分は正常に計算出来ている
・途中から凹凸が出始め,それが成長して発散している
という事態になり,数値的不安定の可能性が強いと思ったからです.オイラー法やルンゲ-クッタ法ではこのような異常振動は進み幅を小さくすれば止まる(ということが数学的に証明されている)ようです.
何卒回答の方,よろしくお願い致します.
お礼
ご回答ありがとうございます.