数値計算法I 第10回連立微分方程式 2006年度九州工業大学工学部電気電子工学コース用講義資料講師：趙孟佑 [email protected] 1 高階微分方程式 • m階の常微分方程式 m 2 m1 d y dy d y d y  f (t, y, , 2 ,L m1 ) m dt dt dt dt これを連立微分方程式に直して、ルンゲ・クッタ法で解く。特に2階微分方程式は多くの物理現象で使う。例:LCR回路 2 高階微分方程式 • m階の常微分方程式 m 2 m1 d y dy d y d y  f (t, y, , 2 ,L m1 ) m dt dt dt dt z1 (t)  y(t),z2 (t)  y(t),L ,zm (t)  y(m1) (t) とおくと、 dz1 dz2 dzm1 dzm (m1)  y  z2 ,  y  z3,L , y  zm ,  y(m) dt dt dt dt である。 3 連立微分方程式 dmy dy d 2 y d m1y  f (t, y, , 2 ,L m1 ) m dt dt dt dt を解くことは、下の1階連立微分方程式を解くに等しい  dz1  dt  z2   dz2  z3  dt  M  dz  m1  zm  dt  dzm (m)  y  f (t, z1, z2 , z3,L , zm )  dt 4 連立微分方程式  dy1  dt  f1 (t, y1, y2 ,L , ym )   dy2  f2 (t, y1, y2 ,L , ym )  dt  M  dy  m1  fm1 (t, y1, y2 ,L , ym )  dt  dym  dt  fm (t, y1, y2 ,L , ym ) t=t0での初期条件は y1 (t0 )  b1, y2 (t0 )  b2 ,L , ym1 (t0 )  bm1, ym (t0 )  bm 5 連立微分方程式の例交流電界中に置かれた電子の速度： dv e   E0 sint dt me 電子の座標x(t)は dx v dt よって、t=t0での初期状態をx(t0)=x0,v(t0)=v0として、  dx  dt  v  dv    e E0 sint me  dt を解いてx(t),v(t)を求める。 6 2階連立微分方程式のルンゲ・クッタ法による解法  dy  dt  f (t, y, z)   dz  g(t, y, z)  dt t  ti で y  yi ,z  zi とする。 t=ti+1=ti+hでの値、yi+1,zi+1は次式で与えられる。 h  y  y   i1 i 6 k1  2k2  2k3  k4    zi1  zi  h l1  2l2  2l3  l4   6 7 2階連立微分方程式のルンゲ・クッタ法による解法 k1  f (ti , yi , zi ) l  g(t , y , z ) i i i 1 h h h  k  f (t  , y  k , z  l1 ) i i 1 i  2 2 2 2  l2  g(ti  h , yi  h k1, zi  h l1 )  2 2 2  k3  f (ti  h , yi  h k2 , zi  h l2 ) 2 2 2   h h h l  g(t  , y  k , z  l2 ) 3 i i 2 i 2 2 2  k4  f (ti  h, yi  hk3, zi  hl3 ) l  g(t  h, y  hk , z  hl ) i i 3 i 3 4 8 2階常微分方程式の例 d2 x dx  3  2x  0 2 dt dt 初期条件はt=0で、x=0,dx/dt=-1 1階の連立微分方程式に直す。  dx  dt  y   dy  3y  2x  dt 初期条件はt=0で、x=0,y=-1 刻み幅h=0.2でt=0から1.0まで計算 x(t)  e2t  et 9 2階常微分方程式の例解析解は t 0.00000000 0.20000000 0.40000000 0.60000000 0.80000000 1.00000000 d2 x dx  3  2x  0 2 dt dt x(t)  e2t  et xi(t) 0.00000000 -0.14833333 -0.22088811 -0.24751482 -0.24734177 -0.23246922 解析解誤差 0.00000000 0.00000000 -0.14841071 0.00007737 -0.22099108 0.00010297 -0.24761742 0.00010260 -0.24743245 0.00009068 -0.23254416 0.00007493 10 2階常微分方程式の例 implicit real*8 (a-h,o-z) parameter(nmax=100000) dimension t(nmax),x(nmax),y(nmax) real*8 k1,k2,k3,k4 real*8 l1,l2,l3,l4 c c h=0.2d0 n=1.0d0/h t(1)=0 x(1)=0 y(1)=-1 do i=1,n t1=t(i) x1=x(i) y1=y(i) k1=f(t1,x1,y1) l1=g(t1,x1,y1) c t2=t(i)+h/2 x2=x(i)+h/2*k1 y2=y(i)+h/2*l1 k2=f(t2,x2,y2) l2=g(t2,x2,y2) c t3=t(i)+h/2 x3=x(i)+h/2*k2 y3=y(i)+h/2*l2 k3=f(t3,x3,y3) l3=g(t3,x3,y3) t4=t(i)+h x4=x(i)+h*k3 y4=y(i)+h*l3 k4=f(t4,x4,y4) l4=g(t4,x4,y4) function f(t,x,y) implicit real*8 (a-h,o-z) f=y return end c function g(t,x,y) implicit real*8 (a-h,o-z) g=-2*x-3*y return end c x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6*h y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6*h t(i+1)=t(i)+h end do do i=1,n+1 sol=dexp(-2*t(i))-dexp(-t(i)) write(6,100)t(i),x(i),sol,dabs(x(i)-sol) end do 100 format(4(1x,f15.8)) stop end c 11