¼ø¶È»ñÎÁ/04

¾ïÈùʬÊýÄø¼°¤Îµá²ò †

¤³¤³¤Ç¤Ï¡¢n ³¬¤Î¾ïÈùʬÊýÄø¼°¤Ï°ì³¬¤Î¾ïÈùʬÊýÄø¼°¤ËÊÑ·Á¤·¡¢u(t) = f(u,t) ¤È¤¤¤¦·Á¼°¤È¤·¤Æ¼è¤ê°·¤¦¤³¤È¤Ë¤·¤Æ¤ª¤¯¡¥

¤µ¤Æ¡¢¾ïÈùʬÊýÄø¼°¤Îµá²ò¤Ë¤Ï¡¢Â绨ÇĤ˸À¤Ã¤Æ

  • °ìÃÊË¡: u(t) ¤Î¶á»÷Ãͤò·×»»¤¹¤ë¤¿¤á¤Ë¡¢u(t-¦¤t) ¤Î¥Ç¡¼¥¿¤À¤±¤¢¤ì¤ÐºÑ¤àÊýË¡¡¥
  • ¿ÃÊË¡: u(t) ¤Î¶á»÷Ãͤò·×»»¤¹¤ë¤¿¤á¤Ë¡¢u(t-¦¤t) ¤À¤±¤Ç¤Ï¤Ê¤¯¡¢u(t-2¦¤t), u(t-3¦¤t) ¤Ê¤É¡¢Ê£¿ô¤Î²áµî¥Ç¡¼¥¿¤òɬÍפȤ¹¤ëÊýË¡

¤ÎÆó¼ïÎब¤¢¤ë¡¥

↑

°ìÃÊË¡ †

´Êñ¤À¤¬°ÂÄêÀ­¤ÈÀºÅÙ¤ËÆñ¤¢¤ê¤Î Euler Ë¡¤È¡¢´è·ò¤Ê Runge-Kutta Ë¡¤ÎÆó¤Ä¤ò¼ø¶È¤Ç¤ÏÀâÌÀ¤·¤¿¡¥ °Ê²¼¡¢¥×¥í¥°¥é¥à¤òºÜ¤»¤Æ¤ª¤³¤¦¡¥

↑

¥µ¥ó¥×¥ë¥×¥í¥°¥é¥à: Euler Ë¡ 1 †

¿Í¸ýÌäÂê¤ËÂФ¹¤ë¤«¤Ê¤ê¥·¥ó¥×¥ë¤Ê¥â¥Ç¥ëÊýÄø¼°¤Ç¤¢¤ë du/dt = u(1-u) ¤È¤¤¤¦¾ïÈùʬÊýÄø¼°¤ËÂФ·¤Æ Euler Ë¡¤Ç¤Î¥×¥í¥°¥é¥à¡¥

  1. # Euler method program for the initial value problem (1).

  2. include Math

  3.  

  4. ##### constants #####

  5. # calculation t from 0 to Tmax.

  6. Tmax = 10.0

  7.  

  8. # Delta t

  9. Dt = 1.0/10

  10.  

  11. # Initial value of u

  12. IniU = 0.4

  13.  

  14. ##### functions #####

  15. # Right hand side of ODE

  16. def f(u)

  17.   return u*(1.0-u)

  18. end

  19.  

  20. # update of the approximate solution

  21. def new_by_euler(u)

  22.   return u + Dt * f(u)

  23. end

  24.  

  25. # result

  26. def print_result(u,t)

  27.   printf("%.6f  %.6f\n", t, u)

  28. end

  29.  

  30. ##### main program #####

  31. # initialization

  32. u = IniU

  33. print_result(u, 0.0)

  34.  

  35. # main loop

  36. for i in (1..(Tmax/Dt).ceil )

  37.   u = new_by_euler(u)

  38.   print_result(u, i*Dt)

  39. end
↑

¥µ¥ó¥×¥ë¥×¥í¥°¥é¥à: Euler Ë¡ 2 †

¥µ¥á¤Èµû¤Î¤è¤¦¤Ê¡¢Èï¿©¼Ô-Êá¿©¼Ô¤Î´Ø·¸¤¬¤¢¤ë¤È¤­¤Î¤½¤Î¼ï¤Î¿ô¤ÎÊÑÆ°¤Ë´Ø¤¹¤ë¥·¥ó¥×¥ë¤Ê¥â¥Ç¥ë¾ïÈùʬÊýÄø¼°¤ËÂФ·¤Æ Euler Ë¡¤Ç¤Î¥×¥í¥°¥é¥à¡¥ ÊýÄø¼°¤Ï µûÁêÅö¤Î¼ï¤Îɤ¿ô¤ò u(t), ¥µ¥áÁêÅö¤Î¼ï¤Îɤ¿ô¤ò v(t) ¤È¤¹¤ë¤È
   du/dt = (2-v)u,
   dv/dt = (2u-3)v,
¤È¤¤¤¦Ï¢Î©°ì³¬¾ïÈùʬÊýÄø¼°¤Ç¤¢¤ë¡¥

  1. # Euler method program for the initial value problem (2).

  2. include Math

  3.  

  4. ##### constants #####

  5. # calculation t from 0 to Tmax.

  6. Tmax = 10.0

  7.  

  8. # Delta t

  9. Dt = 1.0/100

  10.  

  11. # Initial value of u

  12. IniU = [4.0, 1.0]

  13.  

  14. ##### functions #####

  15. # Right hand side of ODE

  16. def f(u)

  17.   a = u[0]

  18.   b = u[1]

  19.   return [ (2.0-b)*a, (2.0*a-3.0)*b ]

  20. end

  21.  

  22. # update of the approximate solution

  23. def new_by_euler(u)

  24.   v = f(u)

  25.   return [ u[0] + Dt*v[0], u[1] + Dt*v[1] ]

  26. end

  27.  

  28. # result

  29. def print_result(u,t)

  30.   printf("%.4f  %.4f  %.4f\n", t, u[0], u[1])

  31. end

  32.  

  33. ##### main program #####

  34. # initialization

  35. u = IniU

  36. print_result(u, 0.0)

  37.  

  38. # main loop

  39. for i in (1..(Tmax/Dt).ceil )

  40.   u = new_by_euler(u)

  41.   print_result(u,i*Dt)

  42. end
↑

¥µ¥ó¥×¥ë¥×¥í¥°¥é¥à: Runge-Kutta Ë¡ 1 †

ƱÍÍ¤Ë du/dt = u(1-u) ¤È¤¤¤¦¾ïÈùʬÊýÄø¼°¤ËÂФ·¤Æ Runge-Kutta Ë¡¤Ç¤Î¥×¥í¥°¥é¥à¡¥

  1. # Runge--Kutta method program for the initial value problem (1).

  2. include Math

  3.  

  4. ##### constants #####

  5. # calculation t from 0 to Tmax.

  6. Tmax = 10.0

  7.  

  8. # Delta t

  9. Dt = 1.0/10

  10.  

  11. # Initial value of u

  12. IniU = 0.4

  13.  

  14. ##### functions #####

  15. # Right hand side of ODE

  16. def f(u)

  17.   return u*(1.0-u)

  18. end

  19.  

  20. # update of the approximate solution

  21. def new_by_rk(u)

  22.   f1 = f(u)

  23.   f2 = f(u + 0.5*Dt*f1)

  24.   f3 = f(u + 0.5*Dt*f2)

  25.   f4 = f(u + Dt*f3)

  26.   return u + Dt*(f1/6.0 + f2/3.0 + f3/3.0 + f4/6.0)

  27. end

  28.  

  29. # result

  30. def print_result(u,t)

  31.   printf("%.6f  %.6f\n", t, u)

  32. end

  33.  

  34. ##### main program #####

  35. # initialization

  36. u = IniU

  37. print_result(u, 0.0)

  38.  

  39. # main loop

  40. for i in (1..(Tmax/Dt).ceil )

  41.   u = new_by_rk(u)

  42.   print_result(u, i*Dt)

  43. end
↑

¥µ¥ó¥×¥ë¥×¥í¥°¥é¥à: Runge-Kutta Ë¡ 2 †

ƱÍͤË
   du/dt = (2-v)u,
   dv/dt = (2u-3)v,
¤È¤¤¤¦Ï¢Î©°ì³¬¾ïÈùʬÊýÄø¼°¤ËÂФ¹¤ë Runge-Kutta Ë¡¤Ç¤Î¥×¥í¥°¥é¥à¡¥

  1. # Runge--Kutta method program for the initial value problem (2).

  2. include Math

  3.  

  4. ##### constants #####

  5. # calculation t from 0 to Tmax.

  6. Tmax = 10.0

  7.  

  8. # Delta t

  9. Dt = 1.0/100

  10.  

  11. # Initial value of u

  12. IniU = [4.0, 1.0]

  13.  

  14. ##### functions #####

  15. # Right hand side of ODE

  16. def f(u)

  17.   a = u[0]

  18.   b = u[1]

  19.   return [ (2.0-b)*a, (2.0*a-3.0)*b ]

  20. end

  21.  

  22. # update of the approximate solution

  23. def new_by_rk(u)

  24.   f1 = f(u)

  25.   f2 = f([ u[0]+0.5*Dt*f1[0], u[1]+0.5*Dt*f1[1] ])

  26.   f3 = f([ u[0]+0.5*Dt*f2[0], u[1]+0.5*Dt*f2[1] ])

  27.   f4 = f([ u[0]+Dt*f3[0], u[1]+Dt*f3[1] ])

  28.   a = u[0] + Dt * (f1[0]/6.0 + f2[0]/3.0 + f3[0]/3.0 + f4[0]/6.0)

  29.   b = u[1] + Dt * (f1[1]/6.0 + f2[1]/3.0 + f3[1]/3.0 + f4[1]/6.0)

  30.   return [ a,b ]

  31. end

  32.  

  33. # result

  34. def print_result(u,t)

  35.   printf("%.4f  %.4f  %.4f\n", t, u[0], u[1])

  36. end

  37.  

  38. ##### main program #####

  39. # initialization

  40. u = IniU

  41. print_result(u, 0.0)

  42.  

  43. # main loop

  44. for i in (1..(Tmax/Dt).ceil )

  45.   u = new_by_rk(u)

  46.   print_result(u,i*Dt)

  47. end
↑

¿ÃÊË¡ †

°ìÃÊË¡¤Î Runge-Kutta Ë¡¤ÈƱÍͤËÀºÅÙ¤ò(ÍýÏÀŪ¤Ë¤Ï)¼«Í³¤Ë¤¢¤²¤é¤ì¤ë¿ÃÊË¡¤È¤·¤Æ¡¢Àþ·Á¿Ãʳ¬Ë¡¤ò¼ø¶È¤ÇÀâÌÀ¤·¤¿¡¥ Àþ·Á¿Ãʳ¬Ë¡¤Ï¡¢·×»»ÀºÅÙ¤ò¤¢¤²¤Æ¤âÈ¡¿ô·×»»¤È¤·¤Æ¤Î·×»»Î̤¬Áý¤¨¤Ê¤¤¤È¤¤¤¦ÂçÊѤ¹¤°¤ì¤¿À­¼Á¤¬¤¢¤ë¡¥ ¤¿¤À¤·¡¢²ò¤ÎÈ¡¿ô¤¬¤¢¤ëÄøÅÙ³ê¤é¤«¤Ç¤¢¤ë¤³¤È¤¬Á°Äó¤È¤µ¤ì¤ë¤³¤È¤ä¡¢²áµî¤Î°ìÄê¤Î¥Ç¡¼¥¿¤ÎÊݸ¤¬É¬Íפʤ³¤È¡¢»þ´ÖÎ¥»¶Éý ¦¤t ¤ò¤½¤Î¾ì¤½¤Î¾ì¤ÇÊѹ¹¤·¤è¤¦¤È¤¹¤ë¤ÈÂçÊѤʤ³¤È¡¢¤È¤¤¤¦¤¢¤¿¤ê¤¬Âå¤ï¤ê¤ËÍ׵ᤵ¤ì¤ë¡¥

¤³¤ì¤Ë¤Ä¤¤¤Æ¤â¥×¥í¥°¥é¥àÎã¤ò¤¢¤²¤Æ¤ª¤³¤¦¡¥

↑

¥µ¥ó¥×¥ë¥×¥í¥°¥é¥à 1 †

¾å¤ÈƱÍÍ¤Ë du/dt = u(1-u) ¤È¤¤¤¦¾ïÈùʬÊýÄø¼°¤ËÂФ¹¤ëÀþ·Á¿Ãʳ¬Ë¡¤Ç¤Î¥×¥í¥°¥é¥à¡¥

  1. # Linear multistep method with 4 stages to the initial value problem 1).

  2. include Math

  3.  

  4. ##### constants #####

  5. Tmax = 10.0

  6. Dt = 1.0/10

  7. IniU = 0.4

  8.  

  9. ##### functions #####

  10. # Right hand side of ODE

  11. def f(u)

  12.   return u*(1.0-u)

  13. end

  14.  

  15. # update of the approximate solution (Runge-Kutta)

  16. def new_by_rk(u)

  17.   f1 = f(u)

  18.   f2 = f(u + 0.5*Dt*f1)

  19.   f3 = f(u + 0.5*Dt*f2)

  20.   f4 = f(u + Dt*f3)

  21.   return u + Dt*(f1/6.0 + f2/3.0 + f3/3.0 + f4/6.0)

  22. end

  23.  

  24. # predictor

  25. def predictor(old_u,old_f)

  26.   return old_u + Dt * ((55.0/24.0)*old_f[0] - (59.0/24.0)*old_f[1] + (37.0/24.0)*old_f[2] - (3.0/8.0)*old_f[3])

  27. end

  28.  

  29. # corrector

  30. def corrector(old_u,old_f,new_f)

  31.   return old_u + Dt * ((3.0/8.0)*new_f + (19.0/24.0)*old_f[0] - (5.0/24.0)*old_f[1] + (1.0/24.0)*old_f[2])

  32. end

  33.  

  34. # result

  35. def print_result(u,t)

  36.   printf("%.6f  %.6f\n", t, u)

  37. end

  38.  

  39. ##### main program #####

  40. # initialization

  41. u = IniU

  42. old_f = [f(u)]

  43. print_result(u,0.0)

  44.  

  45. # we make more three initial values by Runge-Kutta method

  46. for i in (1..3)

  47.   u = new_by_rk(u)

  48.   old_f.unshift(f(u))

  49.   print_result(u,i*Dt)

  50. end

  51.  

  52. # main loop by a LM with PECE cycle

  53. for i in (4..N = (Tmax/Dt).ceil)

  54.   new_u = predictor(u,old_f)       # P: prediction

  55.   new_f = f(new_u)                 # E: Evaluation

  56.   new_u = corrector(u,old_f,new_f) # C: Correction

  57.   new_f = f(new_u)                 # E: Evaluation

  58.  

  59.   u = new_u

  60.   old_f.unshift(new_f)

  61.   old_f.delete_at(4)

  62.   print_result(u,i*Dt)

  63. end
↑

¥µ¥ó¥×¥ë¥×¥í¥°¥é¥à 2 †

ƱÍͤË
   du/dt = (2-v)u,
   dv/dt = (2u-3)v,
¤È¤¤¤¦Ï¢Î©°ì³¬¾ïÈùʬÊýÄø¼°¤ËÂФ¹¤ëÀþ·Á¿Ãʳ¬Ë¡¤Ç¤Î¥×¥í¥°¥é¥à¡¥

  1. # Linear multistep method with 4 stages to the initial value problem 2).

  2. include Math

  3.  

  4. ##### constants #####

  5. Tmax = 10.0

  6. Dt = 1.0/100

  7. IniU = [4.0, 1.0]

  8.  

  9. ##### functions #####

  10. # Right hand side of ODE

  11. def f(u)

  12.   a = u[0]

  13.   b = u[1]

  14.   return [ (2.0-b)*a, (2.0*a-3.0)*b ]

  15. end

  16.  

  17. # update of the approximate solution (Runge-Kutta)

  18. def new_by_rk(u)

  19.   f1 = f(u)

  20.   f2 = f([ u[0]+0.5*Dt*f1[0], u[1]+0.5*Dt*f1[1] ])

  21.   f3 = f([ u[0]+0.5*Dt*f2[0], u[1]+0.5*Dt*f2[1] ])

  22.   f4 = f([ u[0]+Dt*f3[0], u[1]+Dt*f3[1] ])

  23.   a = u[0] + Dt * (f1[0]/6.0 + f2[0]/3.0 + f3[0]/3.0 + f4[0]/6.0)

  24.   b = u[1] + Dt * (f1[1]/6.0 + f2[1]/3.0 + f3[1]/3.0 + f4[1]/6.0)

  25.   return [ a,b ]

  26. end

  27.  

  28. # predictor

  29. def predictor(old_u,old_f)

  30.   a = old_u[0] + Dt * ((55.0/24.0)*old_f[0][0] - (59.0/24.0)*old_f[1][0] + (37.0/24.0)*old_f[2][0] - (3.0/8.0)*old_f[3][0])

  31.   b = old_u[1] + Dt * ((55.0/24.0)*old_f[0][1] - (59.0/24.0)*old_f[1][1] + (37.0/24.0)*old_f[2][1] - (3.0/8.0)*old_f[3][1])

  32.   return [a,b]

  33. end

  34.  

  35. # corrector

  36. def corrector(old_u,old_f,new_f)

  37.   a = old_u[0] + Dt * ((3.0/8.0)*new_f[0] + (19.0/24.0)*old_f[0][0] - (5.0/24.0)*old_f[1][0] + (1.0/24.0)*old_f[2][0])

  38.   b = old_u[1] + Dt * ((3.0/8.0)*new_f[1] + (19.0/24.0)*old_f[0][1] - (5.0/24.0)*old_f[1][1] + (1.0/24.0)*old_f[2][1])

  39.   return [a,b]

  40. end

  41.  

  42. # result

  43. def print_result(u,t)

  44.   printf("%.4f  %.4f  %.4f\n", t, u[0], u[1])

  45. end

  46.  

  47. ##### main program #####

  48. # initialization

  49. u = IniU

  50. old_f = [f(u)]

  51. print_result(u,0.0)

  52.  

  53. # we make more three initial values by Runge-Kutta method

  54. for i in (1..3)

  55.   u = new_by_rk(u)

  56.   old_f.unshift(f(u))

  57.   print_result(u,i*Dt)

  58. end

  59.  

  60. # main loop by a LM with PECE cycle

  61. for i in (4..(Tmax/Dt).ceil)

  62.   new_u = predictor(u,old_f)       # P: prediction

  63.   new_f = f(new_u)                 # E: Evaluation

  64.   new_u = corrector(u,old_f,new_f) # C: Correction

  65.   new_f = f(new_u)                 # E: Evaluation

  66.  

  67.   u = new_u

  68.   old_f.unshift(new_f)

  69.   old_f.delete_at(4)

  70.   print_result(u,i*Dt)

  71. end