Project 1.3_solutions
system:sage


<p>&nbsp;</p>
<div id="_mcePaste" style="position: absolute; left: -10000px; top: 0px; width: 1px; height: 1px; overflow: hidden;">function [f,gradf]=objfun(x)</div>
<div id="_mcePaste" style="position: absolute; left: -10000px; top: 0px; width: 1px; height: 1px; overflow: hidden;">f=(x(1)^2+x(2)^2)^2-x(1)^2-x(2)+x(3)^2;&nbsp;</div>
<div id="_mcePaste" style="position: absolute; left: -10000px; top: 0px; width: 1px; height: 1px; overflow: hidden;">gradf=[4*x(1)*(x(1)^2+x(2)^2)-2*x(1);4*x(2)*(x(1)^2+x(2)^2)-1;2*x(3)];</div>
<div id="_mcePaste" style="position: absolute; left: -10000px; top: 0px; width: 1px; height: 1px; overflow: hidden;">options=optimset('GradObj','on');&nbsp;</div>
<div id="_mcePaste" style="position: absolute; left: -10000px; top: 0px; width: 1px; height: 1px; overflow: hidden;">x0=[1;1;1];[x,fval]=fminunc('objfun',x0,options)</div>
<p><span style="color: #ff0000;">function [f,gradf]=objfun(x)</span></p>
<p><span style="color: #ff0000;">f=(x(1)^2+x(2)^2)^2-x(1)^2-x(2)+x(3)^2;&nbsp;</span></p>
<p><span style="color: #ff0000;">gradf=[4*x(1)*(x(1)^2+x(2)^2)-2*x(1);4*x(2)*(x(1)^2+x(2)^2)-1;2*x(3)];</span></p>
<p><span style="color: #ff0000;">options=optimset('GradObj','on');&nbsp;</span></p>
<p><span style="color: #ff0000;">x0=[1;1;1];[x,fval]=fminunc('objfun',x0,options)</span></p>
<p>&nbsp;</p>

<p>Python code for the "red code"</p>

{{{id=6|
import numpy
from scipy.optimize import fmin_bfgs
def f(x):
    return (x[0]**2+x[1]**2)**2-x[0]**2-x[1]+x[2]**2
def fprime(x):
    grad = numpy.zeros((len(x),), float)
    grad[0] = 4*x[0]*(x[0]**2+x[1]**2)-2*x[0]
    grad[1] = 4*x[1]*(x[0]**2+x[1]**2)-1
    grad[2] = 2*x[2]
    return grad
x0=[1,1,1]

print f(x0), fprime(x0)
#fmin with gradient provided
time xopt = fmin_bfgs(f,x0,fprime) 
print xopt
#fmin with gradient estimated by fmin module
time xopt1 = fmin_bfgs(f,x0) 
print xopt1
///
3 [ 6.  7.  2.]
Optimization terminated successfully.
         Current function value: -0.500000
         Iterations: 11
         Function evaluations: 15
         Gradient evaluations: 15
Time: CPU 0.01 s, Wall: 0.01 s
[ -4.99999944e-01   5.00000405e-01  -3.84579279e-07]
Optimization terminated successfully.
         Current function value: -0.500000
         Iterations: 11
         Function evaluations: 75
         Gradient evaluations: 15
Time: CPU 0.01 s, Wall: 0.01 s
[ -4.99999979e-01   5.00000380e-01  -3.92121138e-07]
}}}

<p>Finds the extreme point analytically and computes the eigenvalues of the Hessian (second derivative in one dimension) to determine</p>
<p>if the function has a minimum or maximum at this point.&nbsp;</p>

{{{id=7|
f(x,y)=x^2*y+y^2+y
f.diff() # gradient
res=solve(list(f.diff()),[x,y])
print res
H=f.diff(2)  # Hessian matrix
print H
H1=H(x=0,y=-1/2)
print H1
H(x=0,y=-1/2).eigenvalues()
///
[
[x == -I, y == 0],
[x == I, y == 0],
[x == 0, y == (-1/2)]
]
[(x, y) |--> 2*y (x, y) |--> 2*x]
[(x, y) |--> 2*x   (x, y) |--> 2]
[-1  0]
[ 0  2]
[-1, 2]
}}}

<p>Sage implementation of the "red code"&nbsp;</p>
<p>1) it computes analytically the extreme point and 2) numerically using the method of BFGS which is a modified Newton</p>

{{{id=8|
sage: var('x , y, z')
sage: f(x,y,z)=(x^2+y^2)^2-x^2-y+z^2 
sage: f.diff() # gradient
sage: res=solve(list(f.diff()),[x,y,z])
sage: print 'roots of gradient', res
sage: x0=[1, 1, 1]
sage: xopt=minimize(f,x0,algorithm="bfgs")
sage: print 'values f at local minimum', f(xopt[0],xopt[1],xopt[2])
sage: H=f.diff(2)  # Hessian matrix
sage: print 'Hessian', H
sage: print 'value of Hessian at xopt', H(xopt[0],xopt[1],xopt[2])
sage: H(xopt[0],xopt[1],xopt[2]).eigenvalues()
///
roots of gradient [
[x == (-2.8551524337e-09 - 9.50618418403e-09*I), y == (-0.314980262474 - 0.545561817986*I), z == 0.0],
[x == (2.8551524337e-09 + 9.50618418403e-09*I), y == (-0.314980262474 - 0.545561817986*I), z == 0.0],
[x == (2.8551524337e-09 - 9.50618418403e-09*I), y == (-0.314980262474 + 0.545561817986*I), z == 0.0],
[x == (-2.8551524337e-09 + 9.50618418403e-09*I), y == (-0.314980262474 + 0.545561817986*I), z == 0.0],
[x == 0.0, y == 0.629960578187, z == 0.0],
[x == (-1/2), y == (1/2), z == 0],
[x == (1/2), y == (1/2), z == 0]
]
Optimization terminated successfully.
         Current function value: -0.500000
         Iterations: 11
         Function evaluations: 15
         Gradient evaluations: 15
values f at local minimum -0.5
Hessian [(x, y, z) |--> 12*x^2 + 4*y^2 - 2              (x, y, z) |--> 8*x*y                  (x, y, z) |--> 0]
[             (x, y, z) |--> 8*x*y     (x, y, z) |--> 4*x^2 + 12*y^2                  (x, y, z) |--> 0]
[                 (x, y, z) |--> 0                  (x, y, z) |--> 0                  (x, y, z) |--> 2]
value of Hessian at xopt [ 2.00000094466 -2.00000139411              0]
[-2.00000139411  4.00000463179              0]
[             0              0              2]
[-13/1383966880035651535063*sqrt(56667688569868626270919004239662009317529) + 345565067138652/115188248656069, 13/1383966880035651535063*sqrt(56667688569868626270919004239662009317529) + 345565067138652/115188248656069, 2]
}}}

<p>3d plotting example</p>

{{{id=10|
var('x,y')
p3d = plot3d(x^2-y^2, (x, -10, 10), (y, -10, 10), mesh=true,frame=true) 
p3d.show()
///
}}}

{{{id=11|
help(minimize)
///
<html><table notracebacks bgcolor="#386074" cellpadding=10 cellspacing=10><tr><td bgcolor="#f5f5f5"><font color="#37546d">
&nbsp;&nbsp;&nbsp;<a target='_new' href='cell://docs-0.html'>Click to open help window</a>&nbsp;&nbsp;&nbsp;
<br></font></tr></td></table></html>
}}}

<p>Implementation of the gradient method for one variable objective function</p>

{{{id=12|
# From calculation, we expect that the local minimum occurs at x=9/4
 
x_old = 0
x_new = 6 # The algorithm starts at x=6
eps = 0.01 # step size
precision = 0.00001
 
def f_prime(x):
    return 4 * x**3 - 9 * x**2
 
while abs(x_new - x_old) > precision:
    x_old = x_new
    x_new = x_old - eps * f_prime(x_old)
print "Local minimum occurs at ", x_new
///
Local minimum occurs at  2.24996460742785
}}}

<p>vector fields of a parametric function</p>

{{{id=13|
var('x,y')
plot_vector_field((sin(x),cos(y)),(x,-pi,pi),(y,-pi,pi))
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=14|
f = e^x + e^y*sin(x)
f.gradient()
///
(e^y*cos(x) + e^x, e^y*sin(x))
}}}

{{{id=15|
g = x^2 + x*y
plot_vector_field(g.gradient(),(x,-2,2),(y,-2,2))
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=16|
var('x,y,z')
G = (x+y,x*y,z)
plot_vector_field3d(G,(x,0,1),(y,0,1),(z,0,1))
///
}}}

{{{id=17|
h = x*y*z
plot_vector_field3d(h.gradient(),(x,-2,2),(y,-2,2),(z,-2,2))
///
}}}

{{{id=18|
var('x,y,z,t')
H = (z,y,x)
r = (cos(6*t),sin(6*t),t)
vfield = plot_vector_field3d(H,(x,-1,1),(y,-1,1),(z,0,pi),colors='orange')
curve = parametric_plot3d(r,(t,0,pi),color='red')
vfield + curve
///
}}}

{{{id=19|
var('theta,phi')
parametric_plot3d((2*sin(phi)*cos(theta),2*sin(phi)*sin(theta),\
2*cos(phi)),(theta,0,2*pi),(phi,0,pi))
///
}}}

{{{id=20|
var('s,t')
parametric_plot3d((s,t,s^2+t^2),(s,-1,1),(t,-1,1))
///
}}}

<p>Gradient method for two variables objective function</p>
<p>Notice that plots the sequence of extreme points on the contour plot of the function</p>

{{{id=22|
var('x,y')
import numpy as np 
xold=vector([1,1])
xnew=vector([0.5, 0.5])
precision=0.00001
pt=[] # defines the list of extreme points it computes in each iteration
eps = 0.01 # step size
f(x,y)=x^2+y^2
fprime=f.gradient()
print(f)
print(fprime)
print "analytic solution of fprime",solve(list(fprime),[x,y])
temp=xnew-xold
nnorm=norm(temp)
plt = contour_plot(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True)
#the implementation of the gradient algorithm
iter=1 
while nnorm > precision:
     xold=xnew
     xnew=vector(xold)-eps*vector(fprime(xold[0],xold[1]))
     nnorm=norm(xnew-xold)
     pt.append(xnew) # sets the list of the sequence of extreme points it computes
     iter=iter+1
print "local minimum at ",xnew
print "number of iterations:",iter
lpt=point(pt)
show(plt+lpt)
///
(x, y) |--> x^2 + y^2
(x, y) |--> (2*x, 2*y)
analytic solution of fprime [
[x == 0, y == 0]
]
local minimum at  (0.000340081752004604, 0.000340081752004604)
number of iterations: 362
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>Computes the exteme points of f(x,y,z) using the method BFGS (a modified Newton)</p>

{{{id=24|
var('x , y, z')
f(x,y,z)=(x^2+y^2)^2-x^2-y+z^2 
fprime=f.diff() # gradient
print(fprime)
res=solve(list(f.diff()),[x,y,z])
print res
x0=[1, 1, 1]
minimize(f,x0,algorithm="bfgs")
///
(x, y, z) |--> (4*(x^2 + y^2)*x - 2*x, 4*(x^2 + y^2)*y - 1, 2*z)
[
[x == (-2.8551524337e-09 - 9.50618418403e-09*I), y == (-0.314980262474 - 0.545561817986*I), z == 0.0],
[x == (2.8551524337e-09 + 9.50618418403e-09*I), y == (-0.314980262474 - 0.545561817986*I), z == 0.0],
[x == (2.8551524337e-09 - 9.50618418403e-09*I), y == (-0.314980262474 + 0.545561817986*I), z == 0.0],
[x == (-2.8551524337e-09 + 9.50618418403e-09*I), y == (-0.314980262474 + 0.545561817986*I), z == 0.0],
[x == 0.0, y == 0.629960578187, z == 0.0],
[x == (-1/2), y == (1/2), z == 0],
[x == (1/2), y == (1/2), z == 0]
]
Optimization terminated successfully.
         Current function value: -0.500000
         Iterations: 11
         Function evaluations: 15
         Gradient evaluations: 15
(-0.499999943818, 0.50000040471, -3.84579279067e-07)
}}}

<p>Gradient method for f(x,y,z)</p>

{{{id=26|
var('x,y,z')
xold=vector([0.5, 0.5, 0.5])
xnew=vector([0.2, 0.2, 0.2])
precision=0.00000001
eps = 0.001 # step size
f(x,y,z)=x^2+y^2+z^2
fprime=f.gradient()
print(f)
print(fprime)
print "analytic solution of fprime",solve(list(fprime),[x,y,z])
temp=xnew-xold
nnorm=norm(temp)
while nnorm > precision:
     xold=xnew
     xnew=vector(xold)-eps*vector(fprime(xold[0],xold[1],xold[2]))
     nnorm=norm(xnew-xold) 
print "local minimum at ",xnew
///
(x, y, z) |--> x^2 + y^2 + z^2
(x, y, z) |--> (2*x, 2*y, 2*z)
analytic solution of fprime [
[x == 0, y == 0, z == 0]
]
local minimum at  (2.87749398873009e-6, 2.87749398873009e-6, 2.87749398873009e-6)
}}}

<p>Computes and plots the extreme point on the contour plot of the function</p>

{{{id=60|

///
}}}

{{{id=29|
var('x, y') 
f = 100 * (y - x^2)^2 + (1 - x)^2 + 100 * (2 - y^2)^2 + (1 - y)^2 
min = minimize(f, [-.1,0.5], disp=1)
print "extreme point",min
plt = contour_plot(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True) 
pt = point(min, color='red', size=50) 
show(plt+pt, aspect_ratio=1, figsize=(4, 4))
pf = plot3d(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True)
show(pf)
///
Optimization terminated successfully.
         Current function value: 4.953484
         Iterations: 8
         Function evaluations: 17
         Gradient evaluations: 17
extreme point (-1.184618747, 1.41254232561)
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=28|
help(minimize)
///
<html><table notracebacks bgcolor="#386074" cellpadding=10 cellspacing=10><tr><td bgcolor="#f5f5f5"><font color="#37546d">
&nbsp;&nbsp;&nbsp;<a target='_new' href='cell://docs-0.html'>Click to open help window</a>&nbsp;&nbsp;&nbsp;
<br></font></tr></td></table></html>
}}}

<p>Displays matrices graphically</p>

{{{id=27|
sage: pi = float(pi)
sage: m = matrix(RDF, 6, [sin(i^2 + j^2) for i in [0,pi/5,..,pi] for j in [0,pi/5,..,pi]])
sage: list_plot3d(m, texture='yellow', frame_aspect_ratio=[1,1,1/3])
///
}}}

{{{id=30|
sage: l=[]
sage: for i in range(6):
...      for j in range(6):
...         l.append((float(i*pi/5),float(j*pi/5),m[i,j]))
sage: list_plot3d(l,texture='yellow')
///
}}}

<p>Displays lists graphically</p>

{{{id=31|
sage: l=[]
sage: for i in range(-5,5):
...    for j in range(-5,5):
...      l.append((normalvariate(0,1),normalvariate(0,1),normalvariate(0,1)))
sage: list_plot3d(l,interpolation_type='nn',texture='yellow',num_points=100)
///
}}}

<p>Applies gradient method and displays the sequence of approximations of the extreme point in the contour graph of&nbsp;</p>
<p>the objective function</p>

{{{id=32|
var('x,y')
import numpy as np 
xold=vector([-.25,.25])
xnew=vector([-0.1, 0.5])
precision=0.00001
pt=[]
eps = 0.1 # step size
f(x,y)=sin(x^2/2-y^2/4+3)*cos(2*x+1-e^y)
fprime=f.gradient()
print(f)
print(fprime)
amin=solve(list(fprime),[x,y])
print "analytic solution of fprime",amin
temp=xnew-xold
nnorm=norm(temp)
plt = contour_plot(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True)
iter=1 
while nnorm > precision:
     xold=xnew
     xnew=vector(xold)-eps*vector(fprime(xold[0],xold[1]))
     nnorm=norm(xnew-xold)
     pt.append(xnew)
     iter=iter+1
print "local minimum at ",xnew
print "number of iterations:",iter
#print pt
lpt=point(pt,color='red')
show(plt+lpt)
pf = plot3d(f, (x, -2, 2), (y, -2, 2),fill=False,cmap='hsv',labels=True)
show(pf)
///
(x, y) |--> sin(1/2*x^2 - 1/4*y^2 + 3)*cos(2*x - e^y + 1)
(x, y) |--> (x*cos(2*x - e^y + 1)*cos(1/2*x^2 - 1/4*y^2 + 3) - 2*sin(2*x - e^y + 1)*sin(1/2*x^2 - 1/4*y^2 + 3), -1/2*y*cos(2*x - e^y + 1)*cos(1/2*x^2 - 1/4*y^2 + 3) + e^y*sin(2*x - e^y + 1)*sin(1/2*x^2 - 1/4*y^2 + 3))
analytic solution of fprime [
x*cos(-2*x + e^y - 1)*cos(-1/2*x^2 + 1/4*y^2 - 3) - 2*sin(-2*x + e^y - 1)*sin(-1/2*x^2 + 1/4*y^2 - 3),
-1/2*y*cos(-2*x + e^y - 1)*cos(-1/2*x^2 + 1/4*y^2 - 3) + e^y*sin(-2*x + e^y - 1)*sin(-1/2*x^2 + 1/4*y^2 - 3)
]
local minimum at  (0.322615850969378, 1.60235597131258)
number of iterations: 145
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>Application of the gradient method for the obective function $(x-y)/(x^2+y^2+1)$</p>
<p>Notice that it displays the sequence of the extreme point approximations on the contour graph of the function&nbsp;</p>

{{{id=35|
var('x,y')
import numpy as np 
xold=vector([0,0])
xnew=vector([-0.1, 0.5])
precision=0.00001
pt=[]
eps = 0.1 # step size
f(x,y)=(x-y)/(x^2+y^2+1) # function to find maximun
fprime=f.gradient()
print(f)
print(fprime)
amin=solve(list(fprime),[x,y])
print "analytic solution of fprime",amin
temp=xnew-xold
nnorm=norm(temp)
plt = contour_plot(f, (x, -3, 3), (y, -3, 3),fill=False,cmap='hsv',labels=True)
iter=1 
while nnorm > precision:
     xold=xnew
     xnew=vector(xold)+eps*vector(fprime(xold[0],xold[1]))
     nnorm=norm(xnew-xold)
     pt.append(xnew)
     iter=iter+1
print "local minimum at ",xnew
print "number of iterations:",iter
#print pt
lpt=point(pt,color='red')
show(plt+lpt)
pf = plot3d(f, (x, -3, 3), (y, -3, 3),fill=False,cmap='hsv',labels=True)
show(pf)
///
(x, y) |--> (x - y)/(x^2 + y^2 + 1)
(x, y) |--> (-2*(x - y)*x/(x^2 + y^2 + 1)^2 + 1/(x^2 + y^2 + 1), -2*(x - y)*y/(x^2 + y^2 + 1)^2 - 1/(x^2 + y^2 + 1))
analytic solution of fprime [
[x == -1/2*sqrt(2), y == 1/2*sqrt(2)],
[x == 1/2*sqrt(2), y == -1/2*sqrt(2)],
[x == -1/2*I*sqrt(2), y == -1/2*I*sqrt(2)],
[x == 1/2*I*sqrt(2), y == 1/2*I*sqrt(2)]
]
local minimum at  (0.707100597445379, -0.706979236153001)
number of iterations: 117
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=55|
y=function('y',x) 
var('x,yy,z')
f=x*y+x^2-y^3*x
ff=f.diff(x).subs({y.diff(x):z,y:yy})
implicit_plot3d(ff,(x,-3,3),(yy,-3,3),(z,-4,4))
///
}}}

{{{id=56|
var('x,y')
f(x,y)=4 * x**3 - 9 * x**2-y
f.gradient()
#lt=plot3d(f,(x,0,2),(y,0,2))
#show(lt)
#plt = contour_plot(f, (x,0,2),(y,0,2),fill=False,cmap='hsv',labels=True)
#show(plt)
fp=plot_vector_field(f.gradient(),(x,0,2),(y,0,4))
fpl=plot(4*x**3-9*x**2,(x,0,2))
show(fpl+fp)
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}