Ȼ/02

Top / Ȼ / 02

f(x)=0 ε

ʬˡ(褬1λ)

򤬴ޤޤƤϤϰϤȾʬ˾Ƥȿˡ ؿϢ³ǡ򤬴ޤޤƤϤϰ(礭Ƥɤ)ΤäƤ OK.

ץץ

  1. include Math
  2.  
  3. def f(x)
  4.   return sin(x)
  5. end
  6.  
  7. em = 1.0*2**(-53)
  8. criteria = 1.2*em # Because sin(x) is almost -(x-PI) around PI.
  9.  
  10. x_l = 0.7*PI
  11. x_r = 1.5*PI
  12.  
  13. while (x_r - x_l > (2*em)*x_l) do
  14.   x_m = (x_r + x_l)/2.0
  15.   if (f(x_m).abs > criteria) then
  16.     if (f(x_l)*f(x_m) < 0) then
  17.       x_r = x_m
  18.     else
  19.       x_l = x_m
  20.     end
  21.   else
  22.     break
  23.   end
  24. end
  25.  
  26. print("Finished. x: ",x_m,", f(x): ",f(x_m),"\n")

ư˽񤭴ơưˡ

f(x) = 0 ˡg(x) := f(x) + x = x Ȥưư⤯֤롥 κݡg(x) ̾ʤ x_new := g(x_old) ȤȿưضŤ

ץץ

  1. include Math
  2.  
  3. def f(x)
  4.   return exp(-x)-x
  5. end
  6.  
  7. def g(x)
  8.   exp(-x)
  9. end
  10.  
  11. # criterion
  12. cri_f = 1.0e-15
  13.  
  14. # initial data
  15. x = 0.5
  16.  
  17. # iteration number
  18. n = 0
  19.  
  20. while (f(x).abs > cri_f) do
  21.   x = g(x)
  22.   n += 1
  23.   print(n, " ",f(x).abs, "\n")
  24. end

Newton ˡ

(ά)

ץץ

  1. include Math
  2.  
  3. Sqt = sqrt(3.0)
  4.  
  5. def f(x,y)
  6.   f_x = x**2+y**2 -1.0
  7.   f_y = Sqt*x-y
  8.   return [f_x,f_y]
  9. end
  10.  
  11. def one_newton(x,y)
  12.   denom = 2.0*(x+Sqt*y)
  13.   u,v = f(x,y)
  14.   new_x = x - ((u+2*y*v)/denom)
  15.   new_y = y - ((Sqt*u-2*x*v)/denom)
  16.   return [new_x,new_y]
  17. end
  18.  
  19. limit_loop = 100
  20. Em = 1.0*2**(-53)
  21. tol_err = 10.0*Em
  22.  
  23. x,y = 0.5,0.5
  24.  
  25. for i in 0..limit_loop do
  26.   f_x,f_y = f(x,y)
  27.   print(i,", ",x,", ",y,", ",f_x,", ",f_y,"\n")
  28.   if ((f_x.abs < tol_err) && (f_y.abs < tol_err)) then
  29.     break
  30.   end
  31.   x,y = one_newton(x,y)
  32. end
  33.  
  34. # output example.
  35.  
  36. 0, 0.5, 0.5, -0.5, 0.366025403784439
  37. 1, 0.549038105676658, 0.950961894323342, 0.205771365940052, 0.0
  38. 2, 0.502189953469069, 0.869818514459077, 0.0087789974610617, 1.11022302462516e-16
  39. 3, 0.500004774982219, 0.866033674296247, 1.91000200775449e-05, 0.0
  40. 4, 0.5000000000228, 0.86602540382393, 9.12012687592778e-11, 0.0
  41. 5, 0.5, 0.866025403784439, -1.11022302462516e-16, 0.0

ۥȥԡˡ

(ά)

ޤʼץ

  1. # Homotopy method to obtain a solution f(x) = x^3-3x+3 = 0.
  2. include Math
  3.  
  4. # error, parameters and initial value
  5. Em = 1.0*2**(-53)
  6. Tol_err = Em*(10.0**2)
  7.  
  8. # We change the homotopy method to the normal Newton method
  9. # when |t-1| < DS or t > 1.
  10. DS = 0.01
  11.  
  12. # Initial value
  13. Ini_x = 3.0
  14.  
  15. # Approximate distance of one movement on the trajectory.
  16. Dl = 0.05
  17.  
  18. # The target function f and its derivative.
  19. def f(x)
  20.   return x**3 -3.0*x +3.0
  21. end
  22.  
  23. def df(x)
  24.   return 3.0*(x**2) -3.0
  25. end
  26.  
  27. # One iteration of Newton method about the function f.
  28. def newton_one_iteration_f(x)
  29.   return x - f(x)/df(x)
  30. end
  31.  
  32. # Newton method to obtain a numerical solution of f(x)=0.
  33. def newton_f(x)
  34.   y = x
  35.   while f(y).abs >Tol_err do
  36.     y = newton_one_iteration_f(y)
  37.   end
  38.   return y
  39. end
  40.  
  41. # Initial value of the function f.
  42. Fw = f(Ini_x)
  43.  
  44. # Approximate movement direction.
  45. def d(x,t)
  46.   gamma = - df(x)/Fw
  47.   dx = - sqrt(Dl/(1.0+(gamma**2)))
  48.   dt = gamma * dx
  49.   return [dx,dt]
  50. end
  51.  
  52. # The function q to be zero on fine movements and its derivative.
  53. def q(dt_tilde,x,t,dx,dt)
  54.   zeta = - dt/dx
  55.   return f(x+dx+zeta*dt_tilde)+(t+dt+dt_tilde - 1.0)*Fw
  56. end
  57.  
  58. def dq(dt_tilde,x,t,dx,dt)
  59.   zeta = - dt/dx
  60.   return df(x+dx+zeta*dt_tilde)*zeta + Fw
  61. end
  62.  
  63. # One iteration of Newton method about the function q.
  64. def newton_one_iteration_q(dt_tilde,x,t,dx,dt)
  65.   return dt_tilde - q(dt_tilde,x,t,dx,dt)/dq(dt_tilde,x,t,dx,dt)
  66. end
  67.  
  68. # Newton method to obtain a numerical solution of q(...)=0.
  69. def newton_q(dt_tilde,x,t,dx,dt)
  70.   y = dt_tilde
  71.   while q(y,x,t,dx,dt).abs >Tol_err do
  72.     y = newton_one_iteration_q(y,x,t,dx,dt)
  73.   end
  74.   return y
  75. end
  76.  
  77. ###################################
  78. # Via a fixed homotopy method we compute numerical solutions
  79. # "x" from t = 0 to t = 1. The last obtained "x" is a numerical
  80. # solution of f(x)=0.
  81.  
  82. t = 0.0
  83. x = Ini_x
  84. print(t," ",x,"\n")
  85.  
  86. while (1.0-t) > DS do
  87.   dx,dt = d(x,t)
  88.   zeta = - dt/dx
  89.   dt_tilde = newton_q(0.0, x,t,dx,dt)
  90.   dx_tilde = zeta*dt_tilde
  91.   x += dx + dx_tilde
  92.   t += dt + dt_tilde
  93.   print(t," ",x,"\n")
  94. end
  95.  
  96. # We use normal Newton method when (1.0-t) <= DS.
  97. x = newton_f(x)
  98. print "1.00000000000000 ",x,"\n"

η

homotopy2-dat.png


źեե: filehomotopy2-dat.png 467 [ܺ]