Comments
Description
Transcript
オイラー法
P83 1階の常微分方程式の初期値問題 1階の常微分方程式 dz = f ( x, z ) dx z ′ = f ( x, z ) 初期条件( x=0のときz=a)を解くには, z ( xk +1 ) = z ( xk ) + xk +1 ∫x f ( x, z )dx k k = 0,1, 2, xj=x0+j h=xj-1+h に zj= z(xj) (ただし,j=0,1,…)を対応させれば z1=z0+h f(x0,z0) z2=z1+h f(x1,z1) : zj+1=zj+h f(xj,zj) これは「微分を差分で置き換えたもの」 で,幾何学的には,「解曲線を,その接 線で置き換えたもの」である.したがっ て,かなりの誤差が生じることになる. j 0 1 2 3 4 5 6 7 8 9 10 x 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 z (オイラー法) 2.00000 1.60000 1.31840 1.10982 0.94970 0.82343 0.72172 0.63838 0.56910 0.51080 0.46123 z (理論値) 2.00000 1.65289 1.38889 1.18343 1.02041 0.88889 0.78125 0.69204 0.61728 0.55402 0.50000 誤差 0.00000 0.05289 0.07049 0.07361 0.07071 0.06546 0.05953 0.05366 0.04818 0.04321 0.03877 相対誤差 0.0% 3.2% 5.1% 6.2% 6.9% 7.4% 7.6% 7.8% 7.8% 7.8% 7.8% 理論解は 2 z 計算 理論 1 0.5 0 0 0.5 1 1.5 x 2 2.5 z(x)のテーラー展開は z(x+h)=z(x)+hz'(x)+z"(x)+… であるが,右辺第3項以下を無視すると z(x+h)≒z(x)+hz'(x) 題意よりz'(x)= f(x,z)だから z(x+h)=z(x)+h f(x,z) となる. 例題2 初期値問題 dz =-xz2 dx ,初期条件(x=1のとき,z=2) を,区間1≦x≦2を10等分して,オイラー法で解け. [解] h=1/10=0.1 x0=1のとき,z0=2 (初期条件) x1=1.1のとき,z1=z0+h f(x0,z0)=2+0.1×{-1.0×22}=1.6 x2=1.2のとき,z2=z1+h f(x1,z1)=1.6+0.1×{-1.1×1.62}=1.3184 x3=1.3のとき,z3=z2+h f(x2,z2) = 1.3184+0.1×{-1.2×1.31842}≒1.10982 x4=1.4のとき,z4=z3+h f(x3,z3) = 1.10982+0.1×{-1.3×1.109822}≒0.94970 : : ルンゲ・クッタ法 オイラー法の改良版 2.5 1.5 オイラー法 dz = − xz 2だから変数分離して dx 1 dz = − xdx を積分して z2 1 2 z −2+1 = − x 2 + c1 z = 2 (∵ c = −2c1 ) −2 + 1 2 x +c 2 2 ∴c = −1 = 0 x = 1のときz = 2だから2 = 1+ c 2 2 ∴z = 2 x z1=z0+hにおけるzの値z1を z1=z0+k ここにk=(k1+2k2+2k3+k4)/6 k1=hf(x0,z0), k2=hf(x0+h/2,z0+k1/2), k3=hf(x0+h/2,z0+k2/2), k4=hf(x0+h,z0+k3) この方法は,解曲線を(x0,z0)のまわりでhに関 して展開したとき,hの4次の項まで,真値と一 致するように公式をきめたことにあたる. 1 2.4.2 連立1階常微分方程式 du1 = f1 (u1 , u2 , u3 , un , t ) dt du2 = f 2 (u1 , u2 , u3 , un , t ) dt : j 0 1 2 3 4 5 6 7 8 9 10 x 1.00000 1.10000 1.20000 1.30000 1.40000 1.50000 1.60000 1.70000 1.80000 1.90000 2.00000 z (オイラー法) 2.00000 1.60000 1.31840 1.10982 0.94970 0.82343 0.72172 0.63838 0.56910 0.51080 0.46123 z理論値 2.00000 1.65289 1.38889 1.18343 1.02041 0.88889 0.78125 0.69204 0.61728 0.55402 0.50000 誤差 0.00000 0.05289 0.07049 0.07361 0.07071 0.06546 0.05953 0.05366 0.04818 0.04321 0.03877 相対誤差 ルンゲクッタ法 誤差 0.0% 2.00000 0.00000 3.2% 1.65292 0.00003 5.1% 1.38892 0.00003 6.2% 1.18346 0.00003 6.9% 1.02044 0.00003 7.4% 0.88891 0.00002 7.6% 0.78127 0.00002 7.8% 0.69206 0.00002 7.8% 0.61730 0.00002 7.8% 0.55403 0.00001 7.8% 0.50001 0.00001 相対誤差 0.000% 0.002% 0.002% 0.003% 0.003% 0.003% 0.003% 0.003% 0.003% 0.002% 0.002% dun = f n (u1 , u2 , u3 dt , un , t ) 初期条件 t=t0のときti,=ai(i=1,2,…, n)の場合に, t (i=1,2,…, n) ui (t k +1 ) = ui (t k ) + ∫t k +1 f i (u1 , u2 , u3 k で近似して,k=0,1,2,……の順に計算 共振(共鳴)現象の微分方程式 2.4.3 高階の常微分方程式 高階(階数をnとする)の常微分方程式は,未知関数 の0階からn-1階までの導関数をそれぞれ別の変数 と置くことにより,n元連立1階常微分方程式 の形に 書き直すことができる.例えば 3 d 2u du +4 + 2u = sin t dt dt 2 du dt , un , t )dt d 2u + 4u = sin 2t dt 2 の一般解は u = c1 cos 2t + c2 sin 2t − t cos 2t 4 ⎧ du ⎪⎪ d t = v ⎨ ⎪ d v = − 4 u + sin 2 t ⎪⎩ d t である. これを数値解法によって解いてみよう。 の場合, = vと置けば du =v dt dv 1 = (sin t − 4v − 2u ) dt 3 という連立1階常微分方程式になる. Sub RungeKutta2() t = 0: 'tの初期値 u = 1: 'uの初期値 v = 1: 'vの初期値 h = 0.1 : 'tの1ステップの変化量 N = 500 : '何回計算するか Cells(3, 1) = "tの初期値" Cells(3, 2) = "uの初期値" Cells(3, 3) = "vの初期値" Cells(3, 4) = "Δt" Cells(3, 5) = "tの終了" Cells(4, 1) = t Cells(4, 2) = u Cells(4, 3) = v Cells(4, 4) = h Cells(4, 5) = t + h * N Cells(6, 2) = "t" Cells(6, 3) = "u" Cells(6, 4) = "v" For i = 0 To N Cells(i + 7, 2).Value = t Cells(i + 7, 3).Value = u Cells(i + 7, 4).Value = v k1 = h * G(t, u, v) l1 = h * F(t, u, v) k2 = h * G(t + h / 2, u + k1 / 2, v + l1 / 2) l2 = h * F(t + h / 2, u + k1 / 2, v + l1 / 2) k3 = h * G(t + h / 2, u + k2 / 2, v + l2 / 2) l3 = h * F(t + h / 2, u + k2 / 2, v + l2 / 2) k4 = h * G(t + h, u + k3, v + l3) l4 = h * F(t + h, u + k3, v + l3) t=t+h u = u + (k1 + 2 * (k2 + k3) + k4) / 6 v = v + (l1 + 2 * (l2 + l3) + l4) / 6 Next i End Sub Function F(t, u, v) 'f=dv/dt F = Sin(2 * t) -4 * u End Function Sub RungeKutta2() t = 0: 'tの初期値 u = 1: 'uの初期値 v = 1: 'vの初期値 h = 0.1 : 'tの1ステップの変化量 N = 500 : '何回計算するか Cells(3, 1) = "tの初期値" Cells(3, 2) = "uの初期値" Cells(3, 3) = "vの初期値" Cells(3, 4) = "Δt" Cells(3, 5) = "tの終了" Cells(4, 1) = t Cells(4, 2) = u Cells(4, 3) = v Cells(4, 4) = h Cells(4, 5) = t + h * N Cells(6, 2) = "t" Cells(6, 3) = "u" Cells(6, 4) = "v" For i = 0 To N Cells(i + 7, 2).Value = t Cells(i + 7, 3).Value = u Cells(i + 7, 4).Value = v k1 = h * G(t, u, v) l1 = h * F(t, u, v) k2 = h * G(t + h / 2, u + k1 / 2, v + l1 / 2) l2 = h * F(t + h / 2, u + k1 / 2, v + l1 / 2) k3 = h * G(t + h / 2, u + k2 / 2, v + l2 / 2) l3 = h * F(t + h / 2, u + k2 / 2, v + l2 / 2) k4 = h * G(t + h, u + k3, v + l3) l4 = h * F(t + h, u + k3, v + l3) t=t+h u = u + (k1 + 2 * (k2 + k3) + k4) / 6 v = v + (l1 + 2 * (l2 + l3) + l4) / 6 Next i End Sub Function F(t, u, v) 'f=dv/dt F = Sin(2 * t) -4 * u End Function Function G(t, u, v) 'g=du/dt=v G=v End Function Function G(t, u, v) 'g=du/dt=v G=v End Function 2