NewtonsMethod.mws
Newton's Multivariable Optimization Method in 2 Variables
by
Dr. William P. Fox (wfox@fmarion.edu) and Dr. William H. Richardson (wrichardson@fmarion.edu)
Department of Mathematics
Francis Marion University, Florence, SC 29501
This application uses Newton's root finding procedure to find a critical point from an initial point. The algorithm uses Cramer's Rule and Newton's Method (see page 73, Mathematical Modeling by Mark Meerschaert, Academic Press, 1993)
Algorithm overview
Variables:
f = the function to be maximized or minimized
x(n) = approximate coordiante of the root after n iterations.
y(n) = approximate coordiante of the root after n iterations.
N = the number of iterations allowed.
F = second partial of f with respect to x
G = second partial of f with respect to y.
INPUTS:
x(0), y(0), N
PROCESS:
Begin
For n = 1 to N do
q <-- Partial of F with respect to x.
r <-- Partial of F with respect to y.
s <--Partial of G with respect to x.
t <--Partial of G with respect to y.
u <-- F(x(n-1),y(n-1))
v <--G(x(n-1),y(n-1))
D <-- ( qt-rs)
x(n) <-- x(n-1) + (ut-vr)/D
y(n) <-- y(n-1) + (qv-su)/D
END
OUTPUT x(n), y(n)
Both the Hessian and the Eigenvalues are displayed to show the test of the critical point for being a maximum, minimum, saddle, or inconclusive.
Initializations
>
restart:with(linalg):
Warning, the protected names norm and trace have been redefined and unprotected
Newton's Program for Optimization
>
Newtons:=proc(f,n::posint,tol::numeric,x11::numeric,x22::numeric)
>
local t1, df1, df2, f1, f2, x, y, dq, dr, ds, dt, A,
q, r, s, t, u, v, newx, newy, xx, yy, dv;
>
t1:=tol;
>
#print(tol);
>
df1:=diff(f,x1);
>
df2:=diff(f,x2);
>
f1:=unapply(df1,x1,x2);
>
f2:=unapply(df2,x1,x2);
>
x:=x11;
>
y:=x22;
>
label_6:
>
dq:=D[1](f1):q:=dq(x,y):#print(q);
>
dr:=D[2](f1):r:=dr(x,y):#print(r);
>
ds:=D[1](f2):s:=ds(x,y):#print(s);
>
dt:=D[2](f2):t:=dt(x,y):#print(t);
>
printf("\nHessian: [ %8.3f %8.3f ]\n",q,r);
>
printf(" [ %8.3f %8.3f ]\n",s,t);
>
A := array([[q,r],[s,t]]);
>
printf("eigenvalues:%8.3f %8.3f\n", evalf(Eigenvals(A))[1],evalf(Eigenvals(A))[2]);
>
A := matrix(2,2, [q,r,s,t]);
>
printf("pos def: %s\n",definite(A, 'positive_def'));
>
u:=-1*f1(x,y):#print(u);
>
v:=-1*f2(x,y):#print(v);
>
dv:=q*t-r*s:#print(dv);
>
newx:=x+((u*t-v*r)/dv);
>
newy:=y+((q*v-s*u)/dv);
>
xx:=x;
>
yy:=y;
>
x:=newx;
>
y:=newy;
>
printf("new x=%8.3f new y=%8.3f\n",newx,newy);
>
if (evalf(sqrt((newx-xx)^2+(newy-yy)^2)) < t1) then
>
printf("\n\nfinal new x=%8.3f final new y=%8.3f\n",newx,newy);
>
printf("final fvalue is %8.3f",evalf(subs({x1=newx,x2=newy},f)));
>
else
>
goto(label_6);
>
fi;
>
#print(newx,newy);
>
end:
>
Examples
>
f:=-(x1-2)^4-(x1-2*x2)^2;
>
Newtons(f,100,.05,0,3);
Hessian: [ -50.000 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -50.378 -7.622
pos def: false
new x= .667 new y= .333
Hessian: [ -23.333 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -24.314 -7.019
pos def: false
new x= 1.111 new y= .556
Hessian: [ -11.481 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -14.103 -5.378
pos def: false
new x= 1.407 new y= .704
Hessian: [ -6.214 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -11.205 -3.009
pos def: false
new x= 1.605 new y= .802
Hessian: [ -3.873 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -10.437 -1.436
pos def: false
new x= 1.737 new y= .868
Hessian: [ -2.832 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -10.178 -.654
pos def: false
new x= 1.824 new y= .912
Hessian: [ -2.370 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -10.076 -.294
pos def: false
new x= 1.883 new y= .941
Hessian: [ -2.164 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -10.033 -.131
pos def: false
new x= 1.922 new y= .961
final new x= 1.922 final new y= .961
final fvalue is -.000
Maximum found at (x,y) =(1.922, 0.961).
>
f:=12*x2-4*x1^2-4*x2^2+4*x1*x2;
>
Newtons(f,100,.01,0.5,1);#4: init (0.5,1)
Hessian: [ -8.000 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -12.000 -4.000
pos def: false
new x= 1.000 new y= 2.000
Hessian: [ -8.000 4.000 ]
[ 4.000 -8.000 ]
eigenvalues: -12.000 -4.000
pos def: false
new x= 1.000 new y= 2.000
final new x= 1.000 final new y= 2.000
final fvalue is 12.000
Maximum found at (x,y) = (1,2).
