授業資料/02 の変更点


* 非線形方程式 f(x)=0 の求解 [#o1cf5870]

** 二分法(定義域が1次元の時) [#dbbbb53f]
解が含まれているはずの範囲を半分に小さくしていく反復法.
関数が連続で、解が含まれているはずの範囲(大きくても良い)を一つ知っていれば OK.

*** サンプルプログラム [#xc3128d8]
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
include Math

def f(x)
  return sin(x)
end

em = 1.0*2**(-53)
criteria = 1.2*em # Because sin(x) is almost -(x-PI) around PI.

x_l = 0.7*PI
x_r = 1.5*PI

while (x_r - x_l > (2*em)*x_l) do
  x_m = (x_r + x_l)/2.0
  if (f(x_m).abs > criteria) then
    if (f(x_l)*f(x_m) < 0) then
      x_r = x_m
    else
      x_l = x_m
    end
  else
    break
  end
end

print("Finished. x: ",x_m,", f(x): ",f(x_m),"\n")
}}


** 不動点問題に書き換えて、不動点を求める方法. [#jd9ce780]

f(x) = 0 を解く代わりに、g(x) := f(x) + x = x という不動点問題の不動点を説く問題に置き換える.
この際、g(x) が縮小写像ならば x_new := g(x_old) という反復で不動点へ近づく.

*** サンプルプログラム [#q59592ff]
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
include Math

def f(x)
  return exp(-x)-x
end

def g(x)
  exp(-x)
end

# criterion
cri_f = 1.0e-15

# initial data
x = 0.5

# iteration number 
n = 0

while (f(x).abs > cri_f) do
  x = g(x)
  n += 1
  print(n, " ",f(x).abs, "\n")
end

}}

** Newton 法 [#lb9d04e5]

(説明略)

*** サンプルプログラム [#p8dd0c72]
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
include Math

Sqt = sqrt(3.0)

def f(x,y)
  f_x = x**2+y**2 -1.0
  f_y = Sqt*x-y
  return [f_x,f_y]
end

def one_newton(x,y)
  denom = 2.0*(x+Sqt*y)
  u,v = f(x,y)
  new_x = x - ((u+2*y*v)/denom)
  new_y = y - ((Sqt*u-2*x*v)/denom)
  return [new_x,new_y]
end

limit_loop = 100
Em = 1.0*2**(-53)
tol_err = 10.0*Em

x,y = 0.5,0.5

for i in 0..limit_loop do
  f_x,f_y = f(x,y)
  print(i,", ",x,", ",y,", ",f_x,", ",f_y,"\n")
  if ((f_x.abs < tol_err) && (f_y.abs < tol_err)) then
    break
  end
  x,y = one_newton(x,y)
end

# output example.

0, 0.5, 0.5, -0.5, 0.366025403784439
1, 0.549038105676658, 0.950961894323342, 0.205771365940052, 0.0
2, 0.502189953469069, 0.869818514459077, 0.0087789974610617, 1.11022302462516e-16
3, 0.500004774982219, 0.866033674296247, 1.91000200775449e-05, 0.0
4, 0.5000000000228, 0.86602540382393, 9.12012687592778e-11, 0.0
5, 0.5, 0.866025403784439, -1.11022302462516e-16, 0.0
}}

** ホモトピー法 [#h077b21c]

(説明略)

*** まじめな実装サンプル [#c6c47ab4]
// programu source 表記
#highlighter(language=ruby,number=on,cache=off){{
# Homotopy method to obtain a solution f(x) = x^3-3x+3 = 0.
include Math

# error, parameters and initial value
Em = 1.0*2**(-53)
Tol_err = Em*(10.0**2)

# We change the homotopy method to the normal Newton method 
# when |t-1| < DS or t > 1.
DS = 0.01

# Initial value
Ini_x = 3.0

# Approximate distance of one movement on the trajectory.
Dl = 0.05

# The target function f and its derivative.
def f(x)
  return x**3 -3.0*x +3.0
end

def df(x)
  return 3.0*(x**2) -3.0
end

# One iteration of Newton method about the function f.
def newton_one_iteration_f(x)
  return x - f(x)/df(x)
end

# Newton method to obtain a numerical solution of f(x)=0.
def newton_f(x)
  y = x
  while f(y).abs >Tol_err do
    y = newton_one_iteration_f(y)
  end
  return y
end

# Initial value of the function f.
Fw = f(Ini_x)

# Approximate movement direction.
def d(x,t)
  gamma = - df(x)/Fw
  dx = - sqrt(Dl/(1.0+(gamma**2)))
  dt = gamma * dx
  return [dx,dt]
end

# The function q to be zero on fine movements and its derivative.
def q(dt_tilde,x,t,dx,dt)
  zeta = - dt/dx
  return f(x+dx+zeta*dt_tilde)+(t+dt+dt_tilde - 1.0)*Fw
end

def dq(dt_tilde,x,t,dx,dt)
  zeta = - dt/dx
  return df(x+dx+zeta*dt_tilde)*zeta + Fw
end

# One iteration of Newton method about the function q.
def newton_one_iteration_q(dt_tilde,x,t,dx,dt)
  return dt_tilde - q(dt_tilde,x,t,dx,dt)/dq(dt_tilde,x,t,dx,dt)
end

# Newton method to obtain a numerical solution of q(...)=0.
def newton_q(dt_tilde,x,t,dx,dt)
  y = dt_tilde
  while q(y,x,t,dx,dt).abs >Tol_err do
    y = newton_one_iteration_q(y,x,t,dx,dt)
  end
  return y
end

###################################
# Via a fixed homotopy method we compute numerical solutions
# "x" from t = 0 to t = 1. The last obtained "x" is a numerical
# solution of f(x)=0.

t = 0.0
x = Ini_x
print(t," ",x,"\n")

while (1.0-t) > DS do
  dx,dt = d(x,t)
  zeta = - dt/dx
  dt_tilde = newton_q(0.0, x,t,dx,dt)
  dx_tilde = zeta*dt_tilde
  x += dx + dx_tilde
  t += dt + dt_tilde
  print(t," ",x,"\n")
end

# We use normal Newton method when (1.0-t) <= DS.
x = newton_f(x)
print "1.00000000000000 ",x,"\n"
}}

*** 上の結果 [#j7b6b74e]

&ref(./homotopy2-dat.png,70%);