Ȼ/04

Top / Ȼ / 04

ʬε

Ǥϡn ξʬϰ쳬ξʬѷu(t) = f(u,t) ȤȤƼ갷ȤˤƤ

ơʬεˤϡ绨Ĥ˸ä

  • ˡ: u(t) ζͤ׻뤿ˡu(t-t) ΥǡкѤˡ
  • ¿ˡ: u(t) ζͤ׻뤿ˡu(t-t) ǤϤʤu(t-2t), u(t-3t) ʤɡʣβǡɬפȤˡ

ब롥

ˡ

ñ٤񤢤 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