Computer Assisted Mathematics (2014)
¼ø¶È»ñÎÁ/02
¤ò¥Æ¥ó¥×¥ì¡¼¥È¤Ë¤·¤ÆºîÀ®
³«»Ï¹Ô:
* ÈóÀþ·ÁÊýÄø¼° f(x)=0 ¤Îµá²ò [#o1cf5870]
** Æóʬˡ(ÄêµÁ°è¤¬1¼¡¸µ¤Î»þ) [#dbbbb53f]
²ò¤¬´Þ¤Þ¤ì¤Æ¤¤¤ë¤Ï¤º¤ÎÈϰϤòȾʬ¤Ë¾®¤µ¤¯¤·¤Æ¤¤¤¯È¿ÉüË¡¡¥
´Ø¿ô¤¬Ï¢Â³¤Ç¡¢²ò¤¬´Þ¤Þ¤ì¤Æ¤¤¤ë¤Ï¤º¤ÎÈÏ°Ï(Â礤¯¤Æ¤âÎɤ¤)...
*** ¥µ¥ó¥×¥ë¥×¥í¥°¥é¥à [#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) arou...
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.20577136594005...
2, 0.502189953469069, 0.869818514459077, 0.00877899746106...
3, 0.500004774982219, 0.866033674296247, 1.91000200775449...
4, 0.5000000000228, 0.86602540382393, 9.12012687592778e-1...
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 meth...
# 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 der...
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...
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 soluti...
# "x" from t = 0 to t = 1. The last obtained "x" is a num...
# 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%);
½ªÎ»¹Ô:
* ÈóÀþ·ÁÊýÄø¼° f(x)=0 ¤Îµá²ò [#o1cf5870]
** Æóʬˡ(ÄêµÁ°è¤¬1¼¡¸µ¤Î»þ) [#dbbbb53f]
²ò¤¬´Þ¤Þ¤ì¤Æ¤¤¤ë¤Ï¤º¤ÎÈϰϤòȾʬ¤Ë¾®¤µ¤¯¤·¤Æ¤¤¤¯È¿ÉüË¡¡¥
´Ø¿ô¤¬Ï¢Â³¤Ç¡¢²ò¤¬´Þ¤Þ¤ì¤Æ¤¤¤ë¤Ï¤º¤ÎÈÏ°Ï(Â礤¯¤Æ¤âÎɤ¤)...
*** ¥µ¥ó¥×¥ë¥×¥í¥°¥é¥à [#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) arou...
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.20577136594005...
2, 0.502189953469069, 0.869818514459077, 0.00877899746106...
3, 0.500004774982219, 0.866033674296247, 1.91000200775449...
4, 0.5000000000228, 0.86602540382393, 9.12012687592778e-1...
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 meth...
# 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 der...
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...
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 soluti...
# "x" from t = 0 to t = 1. The last obtained "x" is a num...
# 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%);
¥Ú¡¼¥¸Ì¾: