* 非線形方程式 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%);
