Chapter 8 M-Files
Script Files
| ShowMySqrt | Plots relative error associated with MySqrt. |
| ShowBisect | Illustrates the method of bisection. |
| ShowNewton | Illustrates the classical Newton iteration. |
| FindConj | Uses fzero to compute Mercury-Earth conjunctions. |
| FindTet | Applies fmin to three different objective functions. |
| ShowGolden | Illustrates Golden section search. |
| ShowSD | Steepest descent test environment. |
| ShowN | Newton test environment. |
| ShowFDN | Finite difference Newton test environment. |
| ShowMixed | Globalized Newton test environment. |
| ShowGN | Gauss-Newton test environment. |
| ShowFmin | Illustrates fmin. |
| ShowFmins | Illustrates fmins. |
Function Files
| MySqrt | Canonical square root finder. |
| Bisection | The method of bisection. |
| StepIsIn | Checks next Newton step. |
| GlobalNewton | Newton method with bisection globalizer. |
| GraphSearch | Interactive graphical search environment. |
| Golden | Golden section search. |
| SDStep | Steepest descent step. |
| NStep | Newton step. |
| FDNStep | Finite difference Newton step. |
| GNStep | Gauss-Newton step. |
| SineMercEarth | The sine of the Mercury-Sun-Earth angle. |
| DistMercEarth | Mercury-Earth distance. |
| TA | Surface area of tetrahedron. |
| TV | Volume of tetrahedron. |
| TE | Total edge length of tetrahedron. |
| TEDiff | Difference between tetrahedron edges. |
| Orbit | Generates and plots orbit points. |
| Sep | Separation between points on two orbits. |
| gSep | Gradient of Sep. |
| gHSep | Gradient and Hessian of sep. |
| SepV | Vector from a point on one orbit to a point on another orbit. |
| JSepV | Jacobian of SepV. |
| rho | Residual of orbit fit. |
| Jrho | Jacobian of rho. |
% Script File: ShowMySqrt
% Plots the error in the function MySqrt.
close all
L = input('Enter left endpoint of test interval:');
R = input('Enter right endpoint of test interval:');
Avals = linspace(L,R,300);
s = zeros(1,300);
for k=1:300
s(k) = MySqrt(Avals(k));
end
error = abs(s-sqrt(Avals))./sqrt(Avals);
figure
plot(Avals,error+eps*.01)
Title('Relative Error in the Function MySqrt(A)')
xlabel('A');
% Script File: ShowBisect
% Illustrates 6 steps of the method of bisection.
close all
fName = input('Enter the name of the function (in quotes):');
a0 = input('Enter the left endpoint of interval of interest:');
b0 = input('Enter the right endpoint of interval of interest:');
% Get initial interval.
x = linspace(a0,b0);
y = feval(fName,x);
figure
plot(x,y,x,0*x);
title('Click in left endpoint a.');
[a,y] = ginput(1);
title('Click in right endpoint b.');
[b,y] = ginput(1);
% Display the initial situation.
x = linspace(a-.1*(b-a),b+.1*(b-a));
y = feval(fName,x);
fa = feval(fName,a);
fb = feval(fName,b);
plot(x,y,x,0*x,[a b], [fa fb],'*')
xlabel(sprintf('At the start: a = %10.6f b = %10.6f',a,b));
for k=1:6
mid = (a+b)/2;
fmid = feval(fName,mid);
title('Strike any key to continue');
pause
if fa*fmid<=0
% There is a root in [a,mid].
b = mid;
fb = fmid;
else
% There is a root in [mid,b].
a = mid;
fa = fmid;
end
x = linspace(a-.1*(b-a),b+.1*(b-a));
y = feval(fName,x);
figure
plot(x,y,x,0*x,[a b], [fa fb],'*')
xlabel(sprintf('After %1.0f steps: a = %10.6f b = %10.6f',k,a,b));
end
% Script File: ShowNewton
% A pure Newton method test environment.
close all
fName = input('Enter the name of the function (with quotes);');
fpName = input('Enter the name of the derivative function (with quotes);');
a = input('Enter left endpoint of interval of interest:');
b = input('Enter right endpoint of interval of interest:');
xvals = linspace(a,b,100);
fvals = feval(fName,xvals);
figure
plot(xvals,fvals,xvals,0*xvals)
v = axis;
title('Click in Starting Value')
[xc,y] = ginput(1);
fc = feval(fName,xc);
fpc = feval(fpName,xc);
Lvals = fc + fpc*(xvals-xc);
hold on
plot(xvals,Lvals,xc,fc,'*')
axis(v)
xlabel(sprintf(' xc = %10.6f fc = %10.3e',xc,fc))
hold off
for k=1:3
step = -fc/fpc;
xc = xc+step;
a = xc-abs(step);
b = xc+abs(step);
xvals = linspace(a,b,100);
fvals = feval(fName,xvals);
fc = feval(fName,xc);
fpc = feval(fpName,xc);
Lvals = fc + fpc*(xvals-xc);
figure
plot(xvals,fvals,xvals,Lvals,xvals,0*xvals,xc,fc,'*')
xlabel(sprintf(' xc = %10.6f fc = %10.3e',xc,fc))
end
% Script File: FindConj
% Estimates spacing between Mercury-Earth Conjunctions
clc
GapEst = input('Spacing Estimate = ');
disp(sprintf('Next ten conjunctions, Spacing Estimate = %8.3f',GapEst))
disp(' ')
t = zeros(11,1);
disp('Conjunction Time Spacing')
disp('-----------------------------------')
for k=1:10;
t(k+1) = fzero('SineMercEarth',k*GapEst);
disp(sprintf(' %2.0f %8.3f %8.3f',k,t(k+1),t(k+1)-t(k)))
end
% Script File: FindTet
% Applies fmin to three different objective functions.
tol = 10e-16;
reps = 20;
disp(' ')
disp('Objective Function Time Theta')
disp('----------------------------------------------------')
t0=clock;
for k=1:reps
OptLat1 = fmin('TA',-pi/2,pi/2,[0 tol]);
end
t1 = etime(clock,t0);
disp(sprintf(' Area %10.3f %20.16f',1,OptLat1))
t0 = clock;
for k=1:reps
OptLat2 = fmin('TV',-pi/2,pi/2,[0 tol]);
end
t2 = etime(clock,t0);
disp(sprintf(' Volume %10.3f %20.16f',t2/t1,OptLat2))
t0 = clock;
for k=1:reps
OptLat3 = fmin('TE',-pi/2,pi/2,[0 tol]);
end
t3 = etime(clock,t0);
disp(sprintf(' Edge %10.3f %20.16f',t3/t1,OptLat3))
OptLat0 = fzero('TEdiff',-.3,eps);
disp(sprintf('\n fzero method: %20.16f',OptLat0))
disp(sprintf(' Exact: %20.16f',asin(-1/3)))
% Script File: ShowGolden
% Illustrates Golden Section Search
close all
tvals = linspace(900,950);
plot(tvals,DistMercEarth(tvals));
axis([900 950 40 150])
axis('off')
hold on
fname = 'DistMercEarth';
a = 900;
b = 950;
r = (3 - sqrt(5))/2;
c = a + r*(b-a); fc = feval(fname,c);
d = a + (1-r)*(b-a); fd = feval(fname,d);
y = 78;
step = 0;
plot([900 950],[y y],[a c d b],[y y y y],'o')
text(952,y,num2str(step))
text(950.5,y+8,'Step')
title('Golden Section Search: (a c d b) Trace')
while step < 5
if fc >= fd
z = c + (1-r)*(b-c);
% [a c d b ] <--- [c d z b]
a = c;
c = d; fc = fd;
d = z; fd = feval(fname,z);
else
z = a + r*(d-a);
% [a c d b ] <--- [a z c d]
b = d;
d = c; fd = fc;
c = z; fc = feval(fname,z);
end
y = y - 7;
step = step+1;
plot([900 950],[y y],[a c d b],[y y y y],'o')
text(952,y,num2str(step))
end
tmin = (c+d)/2;
hold off
% Script File: ShowSD
% Minimizes the sep function using steepest descent.
planet1 = struct('A',10,'P',2,'phi', pi/8);
planet2 = struct('A', 4,'P',1,'phi',-pi/7);
close all
Lmax = 1; % Max line step
nSteps = 25; % Total number of steps.
manualSteps = 3; % Number of manual steps with point-and-click line search.
con_size = 20; % Size of the contour plot
% Plot the two orbits.
figure
t=linspace(0,2*pi);
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');
hold off
% Display contours of the sep function.
tRange = linspace(0,2*pi,con_size);
Z = zeros(con_size,con_size);
for i=1:con_size
for j=1:con_size
t = [tRange(i);tRange(j)];
Z(i,j) = Sep(t,planet1,planet2);
end
end
figure
C = contour(tRange,tRange,Z',10);
title('Enter Contour Labels (Optional). To quit, strike .')
xlabel('t1')
ylabel('t2')
Clabel(C,'manual')
% Get ready for the iteration.
title('Enter Starting Point for Steepest Descent')
tc = zeros(2,1);
[tc(1),tc(2)] = ginput(1);
fc = Sep(tc,planet1,planet2);
gc = gSep(tc,planet1,planet2);
tvals = tc
fvals = fc;
gvals = norm(gc);
title('')
% Steepest Descent with line search:
clc
disp('Step sep t(1) t(2) norm(grad)')
disp('----------------------------------------------------------------')
disp(sprintf('%2.0f %10.6f %10.6f %10.6f %8.3e ',0,fc,tc,norm(gc)))
figure(3)
for step = 1:nSteps
if step<=manualSteps
% Manual Line Search
[tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,0);
else
% Automated Line Search
[tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,1);
end
tvals = [tvals tnew]; tc = tnew;
fvals = [fvals fnew]; fc = fnew;
gvals = [gvals norm(gnew)]; gc = gnew;
disp(sprintf('%2.0f %10.6f %10.6f %10.6f %8.3e ',step,fc,tc,norm(gc)))
end
% Show solution on orbit plot:
figure(1)
hold on
last = nSteps+1;
pt1 = Orbit(tvals(1,last),planet1,'*');
pt2 = Orbit(tvals(2,last),planet2,'*');
plot([pt1.x pt2.x],[pt1.y pt2.y])
title(sprintf('min Sep = %8.4f',fvals(last)))
hold off
% Show the descent path on the contour plot:
figure(2)
hold on
plot(tvals(1,:),tvals(2,:),tvals(1,1),tvals(2,1),'o')
hold off
title(sprintf('tmin = (%8.3f,%8.3f) norm(gmin)= %8.4e',tvals(1,last),tvals(2,last),gvals(last)))
% Plot the descent of the sep and its gradient:
figure(3)
subplot(2,1,1)
plot(fvals)
title('Value of Sep')
xlabel('Iteration')
subplot(2,1,2)
semilogy(gvals)
title('Value of norm(gSep).')
xlabel('Iteration')
% Script File: ShowN
% Newton method test environment.
close all
planet1 = struct('A',15,'P',2,'phi', pi/10);
planet2 = struct('A',20,'P',3,'phi',-pi/8);
itmax = 10;
tol = 1e-14;
% Plot Orbits
figure
axis equal off
t = linspace(0,2*pi);
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');
% Enter Starting Points
tc = rand(2,1)*2*pi;
Orbit(tc(1),planet1,'*');
Orbit(tc(2),planet2,'*');
title('Initial Guess')
hold off
% Initializations
Fc = SepV(tc,planet1,planet2); r = norm(Fc);
Jc = JSepV(tc,planet1,planet2); c = cond(Jc);
clc
disp('Iteration tc(1) tc(2) norm(Fc) cond(Jc)')
disp('------------------------------------------------------------------------------')
disp(sprintf(' %2.0f %20.16f %20.16f %8.3e %8.3e',0,tc(1),tc(2),r,c))
% The Newton Iteration
k=0;
while (ktol)
% Take a step and display
[tc,Fc,Jc] = NStep(tc,Fc,Jc,planet1,planet2);
k=k+1;
figure
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');
Orbit(tc(1),planet1,'*');
Orbit(tc(2),planet2,'*');
hold off
r = norm(Fc);
c = cond(Jc);
title(sprintf(' Iteration = %2.0f norm(Fc) = %8.3e',k,r))
disp(sprintf(' %2.0f %20.16f %20.16f %8.3e %8.3e',k,tc(1),tc(2),r,c))
end
sol = Orbit(tc(1),planet1);
disp(sprintf('\n Intersection = ( %17.16f , %17.16f )',sol.x,sol.y))
% Script File: ShowFDN
% Finite Difference Newton method test environment.
close all
planet1 = struct('A',15,'P',2,'phi', pi/10);
planet2 = struct('A',20,'P',3,'phi',-pi/8);
itmax = 10;
tol = 1e-14;
% Plot Orbits
figure
axis equal off
t = linspace(0,2*pi);
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');
% Enter Starting Points
tc = rand(2,1)*2*pi;
Orbit(tc(1),planet1,'*');
Orbit(tc(2),planet2,'*');
title('Initial Guess')
hold off
% Initializations
Fc = SepV(tc,planet1,planet2); r = norm(Fc);
clc
disp('Iteration tc(1) tc(2) norm(Fc)')
disp('--------------------------------------------------------------------')
disp(sprintf(' %2.0f %20.16f %20.16f %8.3e',0,tc(1),tc(2)))
% The FDNewton Iteration
k=0;
while (ktol)
% Take a step and display
[tc,Fc] = FDNStep(tc,Fc,planet1,planet2);
k=k+1;
figure
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');
Orbit(tc(1),planet1,'*');
Orbit(tc(2),planet2,'*');
hold off
r = norm(Fc);
title(sprintf(' Iteration = %2.0f norm(Fc) = %8.3e',k,r))
disp(sprintf(' %2.0f %20.16f %20.16f %8.3e %8.3e',k,tc(1),tc(2),r))
end
sol = Orbit(tc(1),planet1);
disp(sprintf('\n Intersection = ( %17.16f , %17.16f )',sol.x,sol.y))
% Script File: ShowMixed
% Minimizes the sep function by doing a few steepest descent steps
% (with line search) and then a few full Newton steps.
close all
SDSteps = 2; % The number of steepest descent steps.
NSteps = 10; % The number of Newton steps.
planet1 = struct('A',10,'P',2,'phi', pi/8);
planet2 = struct('A', 4,'P',1,'phi',-pi/7);
Lmax = 1; % Max line step
con_size = 20; % Size of the contour plot
% Plot the two orbits.
figure
t=linspace(0,2*pi);
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');
hold off
% Display contours of the sep function.
tRange = linspace(0,2*pi,con_size);
Z = zeros(con_size,con_size);
for i=1:con_size
for j=1:con_size
t = [tRange(j);tRange(i)];
Z(i,j) = Sep(t,planet1,planet2);
end
end
figure
C = contour(tRange,tRange,Z,10);
title('Enter Contour Labels (Optional). To quit, strike .')
xlabel('t1')
ylabel('t2')
Clabel(C,'manual')
% Get ready for the iteration.
title('Enter Starting Point for Steepest Descent')
tc = zeros(2,1);
[tc(1),tc(2)] = ginput(1);
fc = Sep(tc,planet1,planet2);
gc = gSep(tc,planet1,planet2);
tvals = tc;
fvals = fc;
gvals = norm(gc);
title('')
% Steepest Descent with line search:
clc
disp('Step sep t(1) t(2) norm(grad)')
disp('------------------------------------------------------------------------')
disp(sprintf('%2.0f %10.6f %18.15f %18.15f %5.1e ',0,fc,tc,norm(gc)))
figure(3)
for step = 1:SDSteps
[tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,0);
tvals = [tvals tnew]; tc = tnew;
fvals = [fvals fnew]; fc = fnew;
gvals = [gvals norm(gnew)]; gc = gnew;
disp(sprintf('%2.0f %10.6f %18.15f %18.15f %5.1e ',step,fc,tc,norm(gc)))
end
disp('------------------------------------------------------------------------')
[gc,Hc] = gHSep(tc,planet1,planet2);
for step=SDSteps+1:SDSteps+NSteps;
tc = tc - Hc\gc;
[gc,Hc] = gHSep(tc,planet1,planet2);
fc = Sep(tc,planet1,planet2);
tvals = [tvals tc];
fvals = [fvals fc];
gvals = [gvals norm(gc)];
disp(sprintf('%2.0f %10.6f %18.15f %18.15f %5.1e ',step,fc,tc,norm(gc)))
end
% Show solution on orbit plot:
figure(1)
hold on
last = length(gvals);
pt1 = Orbit(tvals(1,last),planet1,'*');
pt2 = Orbit(tvals(2,last),planet2,'*');
plot([pt1.x pt2.x],[pt1.y pt2.y])
title(sprintf('min Sep = %8.4f',fvals(last)))
hold off
% Show the descent path on the contour plot:
figure(2)
hold on
plot(tvals(1,:),tvals(2,:),tvals(1,1),tvals(2,1),'o')
hold off
title(sprintf('tmin = (%8.3f,%8.3f) norm(gmin)= %8.4e',tvals(1,last),tvals(2,last),gvals(last)))
% Plot the descent of the sep and its gradient:
figure(4)
subplot(2,1,1)
plot(fvals)
title('Value of Sep')
xlabel('Iteration')
subplot(2,1,2)
semilogy(gvals)
title('Value of norm(gSep).')
xlabel('Iteration')
% Script File: ShowGN
% Illustrates Gauss-Newton with line search applied to the problem of fitting a
% two-parameter ellipse to some noisy, nearly-elliptical data.
A0 = 5; P0 = 1; % Parameters of the "true" ellipse.
del = .1; % Noise factor.
m = 50; % Number of points.
con_size = 20; % Size of contour plot
kmax = 25 ; % Maximum number of GN steps.
Lmax = 8; % Maximum step length factor
tol = .000001; % Size of gradient for termination.
% Generate a noisy orbit
theta = linspace(0,2*pi,m)';
c = cos(linspace(0,2*pi,m)');
r = ((2*A0*P0)./(P0*(1-c) + A0*(1+c))).*(1 + del*randn(m,1));
plist = [r c];
% Plot the true orbit and the radius vector samples
close all
figure
axis equal
planet = struct('A',A0,'P',P0,'phi',0);
%Orbit(linspace(0,2*pi,200),planet,'-');
hold on
x = r.*c;
y = r.*sin(theta);
plot(x,y,'o')
hold off
% Generate Contour Plot
range = linspace(.9*min(r),1.5*max(r),con_size);
Z = zeros(con_size,con_size);
for i=1:con_size
for j=1:con_size
f = rho([range(j);range(i)],plist);
Z(i,j) = (f'*f)/2;
end
end
figure
C = contour(range,range,Z,20);
title('Enter Contour Labels (Optional). To quit, strike .')
xlabel('A'); ylabel('P'); Clabel(C,'manual')
% Initialize for iteration
title('Enter Starting Point for Gauss-Newton')
zc = zeros(2,1);
[zc(1),zc(2)] = ginput(1);
Jc = Jrho(zc,plist);
Fc = rho(zc,plist);
wc = (Fc'*Fc)/2;
k=0;
gnorm = norm(Jc'*Fc);
clc
disp(sprintf('"True A" = %5.3f, "True P" = %5.3f',A0,P0))
disp(sprintf('Relative Noise = %5.3f, Samples = %3.0f\n',del,m))
disp('Iteration wc A P norm(grad)')
disp('----------------------------------------------------------')
disp(sprintf(' %2.0f %10.4f %10.6f %10.6f %5.3e ',k,wc,zc(1),zc(2),gnorm))
figure(3)
% Gauss-Newton with linesearch
while (gnorm > tol ) & (k < kmax)
if k<=2
[zc,Fc,wc,Jc] = GNStep('rho','Jrho',zc,Fc,wc,Jc,plist,Lmax,0);
else
[zc,Fc,wc,Jc] = GNStep('rho','Jrho',zc,Fc,wc,Jc,plist,Lmax,1);
end
k=k+1;
gnorm = norm(Jc'*Fc);
disp(sprintf(' %2.0f %10.4f %10.6f %10.6f %5.3e ',k,wc,zc(1),zc(2),gnorm))
end
figure(1)
planetModel = struct('A',zc(1),'P',zc(2),'phi',0);
hold on
Orbit(linspace(0,2*pi),planetModel,'-');
hold off
% Script File: ShowFmin
% Illustrates the 1-d minimizer fmin.
clc
tol = .000001;
disp('Local minima of the Mercury-Earth separation function.')
disp(sprintf('\ntol = %8.3e\n\n',tol))
disp(' Initial Interval tmin f(tmin) f evals')
disp('----------------------------------------------------')
options = zeros(18,1);
for k=1:8
L = 100+(k-1)*112;
R = L+112;
options(2) = tol;
[tmin options] = fmin('DistMercEarth',L,R,options);
minfeval = options(8);
nfevals = options(10);
disp(sprintf(' [%3.0f,%3.0f] %10.5f %10.5f, %6.0f',L,R,tmin,minfeval,nfevals))
end
% Script Files: ShowFmins
% Illustrates the multidimensional minimizer fmins.
close all
planet1 = struct('A',10,'P',2,'phi', pi/8);
planet2 = struct('A', 4,'P',1,'phi',-pi/7);
% Plot the orbots
tVals = linspace(0,2*pi);
Orbit(tVals,planet1,'-');
hold on
Orbit(tVals,planet2,'--');
% Call fmins
trace = 0;
steptol = .000001;
ftol = .000001;
options = [trace steptol ftol];
t0=[5;2];
[t,options] = Fmins('Sep',t0,options,[],planet1,planet2);
% Show Solution
pt1 = Orbit(t(1),planet1,'o');
pt2 = Orbit(t(2),planet2,'o');
plot([pt1.x pt2.x],[pt1.y pt2.y])
hold off
its = options(10);
gval = norm(gSep(t,planet1,planet2));
title(sprintf('t_{*} = ( %8.6f , %8.6f ), Steps = %3.0f, norm(gradient) = %5.3e',t,its,gval))
axis off equal
function x = MySqrt(A)
% x = MySqrt(A)
% A is a non-negative real and x is its square root.
if A==0
x = 0;
else
TwoPower = 1;
m = A;
while m>=1
m = m/4;
TwoPower = 2*TwoPower;
end
while m < .25
m = m*4;
TwoPower = TwoPower/2;
end
% sqrt(A) = sqrt(m)*TwoPower
x = (1+2*m)/3;
for k=1:4
x = (x + (m/x))/2;
end
x = x*TwoPower;
end
function root = Bisection(fname,a,b,delta)
% root = Bisection(fname,a,b,delta)
% Method of bisection.
%
% fname is a string that names a continuous function f(x) of a single variable.
% a and b define an interval [a,b] and f(a)f(b) < 0
% delta non-negative real.
%
% root is the midpoint of an interval [alpha,beta]
% with the property that f(alpha)f(beta)<=0 and
% |
% |beta-alpha| <= delta+eps*max(|alpha|,|beta|)
fa = feval(fname,a);
fb = feval(fname,b);
if fa*fb > 0
disp('Initial interval is not bracketing.')
return
end
if nargin==3
delta = 0;
end
while abs(a-b) > delta+eps*max(abs(a),abs(b))
mid = (a+b)/2;
fmid = feval(fname,mid);
if fa*fmid<=0
% There is a root in [a,mid].
b = mid;
fb = fmid;
else
% There is a root in [mid,b].
a = mid;
fa = fmid;
end
end
root = (a+b)/2;
function ok = StepIsIn(x,fx,fpx,a,b) % ok = StepIsIn(x,fx,fpx,a,b) % Yields 1 if the next Newton iterate is in [a,b] and 0 otherwise. % x is the current iterate, fx is the value of f at x, and fpx is % the value of f' at x. if fpx > 0 ok = ((a-x)*fpx <= -fx) & (-fx <= (b-x)*fpx); elseif fpx < 0 ok = ((a-x)*fpx >= -fx) & (-fx >= (b-x)*fpx); else ok = 0; end
function [x,fx,nEvals,aF,bF] = GlobalNewton(fName,fpName,a,b,tolx,tolf,nEvalsMax)
% [ x,fx,nEvals,aF,bF] = GlobalNewton(fName,fpName,a,b,tolx,tolf,nEvalsMax)
% fName string that names a function f(x).
% fpName string that names the derivative function f'(x).
% a,b A root of f(x) is sought in the interval [a,b]
% and f(a)*f(b)<=0.
% tolx,tolf Nonnegative termination criteria.
% nEvalsMax Maximum number of derivative evaluations.
%
% x An approximate zero of f.
% fx The value of f at x.
% nEvals The number of derivative evaluations required.
% aF,bF The final bracketing interval is [aF,bF].
%
% The iteration terminates as soon as x is within tolx of a true zero or
% if |f(x)|<= tolf or after nEvalMax f-evaluations.
fa = feval(fName,a);
fb = feval(fName,b);
if fa*fb>0
disp('Initial interval not bracketing.')
return
end
x = a;
fx = feval(fName,x);
fpx = feval(fpName,x);
disp(sprintf('%20.15f %20.15f %20.15f',a,x,b))
nEvals = 1;
while (abs(a-b) > tolx ) & (abs(fx) > tolf) & ((nEvals < nEvalsMax) | (nEvals==1))
%[a,b] brackets a root and x = a or x = b.
if StepIsIn(x,fx,fpx,a,b)
%Take Newton Step
disp('Newton')
x = x-fx/fpx;
else
%Take a Bisection Step:
disp('Bisection')
x = (a+b)/2;
end
fx = feval(fName,x);
fpx = feval(fpName,x);
nEvals = nEvals+1;
if fa*fx<=0
% There is a root in [a,x]. Bring in right endpoint.
b = x;
fb = fx;
else
% There is a root in [x,b]. Bring in left endpoint.
a = x;
fa = fx;
end
disp(sprintf('%20.15f %20.15f %20.15f',a,x,b))
end
aF = a;
bF = b;
function [L,R] = GraphSearch(fname,a,b,Save,nfevals)
% [L,R] = GraphSearch(fname,a,b,Save,nfevals)
%
% Graphical search. Produces sequence of plots of the function f(x). The user
% specifies the x-ranges by mouseclicks.
%
% name is a string that names a function f(x) that is defined
% on the interval [a,b].
% nfevals>=2
%
% Save is used to determine how the plots are saved. If Save is
% nonzero, then each plot is saved in a separate figure window.
% If Save is zero or if GraphSearch is called with just three
% arguments, then only the final plot is saved.
%
% [L,R] is the x-range of the final plot. The plots are based on nfevals
% function evaluations. If GraphSearch is called with less
% than five arguments, then nfevals is set to 100.
close all
if nargin==3
Save=0;
end
if nargin<5
nfevals=100;
end
AnotherPlot = 1;
L = a;
R = b;
while AnotherPlot
x = linspace(L,R,nfevals);
y = feval(fname,x);
if Save
figure
end
ymin = min(y);
plot(x,y,[L R],[ymin ymin])
title(['The Function ' fname])
v = axis;
v(1) = L;
v(2) = R;
v(3) = v(3)-(v(4)-v(3))/10;
axis(v)
xlabel('Enter New Left Endpoint. (Click off x-range to terminate.)')
text(R,ymin,[' ' num2str(ymin)])
[x1,y1] = ginput(1);
if (x1 < L) | (R < x1)
AnotherPlot=0;
else
xlabel('Enter New Right Endpoint. (Click off x-range to terminate.)')
[x2,y2] = ginput(1);
if (x2 < L) | (R < x2)
AnotherPlot=0;
else
L = x1;
R = x2;
end
end
xlabel(' ')
end
function tmin = Golden(fname,a,b)
% tmin = Golden(fname,a,b)
% Golden Section Search
%
% fname string that names function f(t) of a single variable.
% a,b define an interval [a,b] upon which f is unimodal.
%
% tmin approximate global minimizer of f on [a,b].
r = (3 - sqrt(5))/2;
c = a + r*(b-a); fc = feval(fname,c);
d = a + (1-r)*(b-a); fd = feval(fname,d);
while (d-c) > sqrt(eps)*max(abs(c),abs(d))
if fc >= fd
z = c + (1-r)*(b-c);
% [a c d b ] <--- [c d z b]
a = c;
c = d; fc = fd;
d = z; fd = feval(fname,z);
else
z = a + r*(d-a);
% [a c d b ] <--- [a z c d]
b = d;
d = c; fd = fc;
c = z; fc = feval(fname,z);
end
end
tmin = (c+d)/2;
function [tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,auto)
% [tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,auto)
% Generates a steepest descent step.
%
% fName the name of a function f(x,plist) that accepts column
% n-vectors and returns a scalar.
% gName the name of a function g(x,plist) that returns the
% gradient of f at x.
% xc a column-vector, an approximate minimizer
% fc f(xc,plist)
% gc g(xc,plist)
% plist orbit parameters for f and g
% Lmax max step length
% auto if auto is nonzero, then automated line search
%
% xnew a column n-vector, an improved minimizer
% fnew f(xnew,plist)
% gnew g(new,plist)
%
nplotvals = 20;
% Get the Steepest Descent Search Direction
sc = -gc;
% Line Search
% Try to get L<=Lmax so xc+L*sc is at least as good as xc.
L = Lmax;
Halvings = 0;
fL = Sep(tc+L*sc,planet1,planet2);
while (fL>=fc) & Halvings<=10
L = L/2;
Halvings = Halvings+1;
fL = Sep(tc+L*sc,planet1,planet2);
end
% Sample across [0.L]
lambdavals = linspace(0,L,nplotvals);
fvals = zeros(1,nplotvals);
for k=1:nplotvals
fvals(k) = Sep(tc+lambdavals(k)*sc,planet1,planet2);
end
if auto==0
% Manual line search.
plot(lambdavals,fvals);
xlabel('lambda');
ylabel('Sep(tc+lambda*sc)');
title('Click the Best lambda value.');
[lambda,y] = ginput(1);
tnew = tc+lambda*sc;
fnew = Sep(tnew,planet1,planet2);
gnew = gSep(tnew,planet1,planet2);
else
% Automated line search.
[fnew,i] = min(fvals(1:nplotvals));
tnew = tc + lambdavals(i(1))*sc;
gnew = gSep(tnew,planet1,planet2);
end
function [tnew,Fnew,Jnew] = NStep(tc,Fc,Jc,planet1,planet2) % [tnew,Fnew,Jnew] = NStep(tc,Fc,Jc,planet1,planet2) % Newton Step % % tc is a column 2-vector, Fc is the value of SepV(t,planet1,planet2) at t=tc % and Jc is the value of JSepV(t,planet1,planet2) at t=tc. Does one Newton step % applied to SepV rendering a new approximate root tnew. Fnew and Jnew are the % values of SepV(tnew,planet1,planet2) and jSepV(tnew,planet1,planet2) respectively. sc = -(Jc\Fc); tnew = tc + sc; Fnew = SepV(tnew,planet1,planet2); Jnew = JSepV(tnew,planet1,planet2);
function [tnew,Fnew] = FDNStep(tc,Fc,planet1,planet2) % [tnew,Fnew] = FDNStep(tc,Fc,planet1,planet2) % Finite difference Newton step. % % tc is a column 2-vector, Fc is the value of SepV(t,planet1,planet2) at t=tc % Does one finite difference Newton step applied to SepV rendering a new % approximate root tnew. Fnew is the values of SepV(tnew,planet1,planet2). % Set up the FD Jacobian. Jc = zeros(2,2); delta = sqrt(eps); for k=1:2 tk = tc; tk(k) = tc(k) + delta; Jc(:,k) = (SepV(tk,planet1,planet2) - Fc)/delta; end % The Finite Difference Newton Step. sc = -(Jc\Fc); tnew = tc + sc; Fnew = SepV(tnew,planet1,planet2);
function [xnew,Fnew,wnew,Jnew] = GNStep(FName,JName,xc,Fc,wc,Jc,plist,Lmax,auto)
% [xnew,Fnew,wnew,Jnew] = GNStep(FName,JName,xc,Fc,wc,Jc,plist,Lmax,auto)
% Generates a Gauss-Newton Step
%
% FName the name of a function F(x,plist) that accepts a column
% n-vectors and column m-vector
% JName the name of a function JF(x,plist) that returns the
% Jacobian of F at x.
% xc a column n-vector, an approximate minimizer.
% Fc F(xc,plist)
% wc Fc'*Fc/2
% Jc JF(xc,plist)
% plist parameter list
% auto 0 for manual line search and nonzero otherwise.
%
% xnew a column n-vector, an improved minimizer.
% Fnew F(xnew,plist)
% wnew Fnew'*Fnew/2
% Jnew JF(xnew,plist)
nplotvals = 20; % Number of F evaluations per line search.
% Get the Gauss-Newton Search Direction
sc = -Jc\Fc;
% Line Search
% Try to get L<=Lmax so xc+L*sc is at least as good as xc.
L = Lmax;
Halvings = 0;
FL = feval(FName,xc+L*sc,plist);
while ((FL'*FL/2)>=wc) & Halvings<=10
L = L/2;
Halvings = Halvings+1;
FL = feval(FName,xc+L*sc,plist);
end
lambdavals = linspace(0,L,nplotvals);
wvals = zeros(nplotvals,1);
for k=1:nplotvals
F = feval(FName,xc+lambdavals(k)*sc,plist);
wvals(k) = (F'*F)/2;
end
if auto==0
% Manual line search
plot(lambdavals,wvals);
xlabel('lambda');
ylabel('w(xc+lambda*sc)');
title('Enter the Best lambda value.');
[lambda,y] = ginput(1);
xnew = xc+lambda*sc;
else
% Automated line search
[fnew,i] = min(wvals(1:nplotvals));
xnew = xc + lambdavals(i(1))*sc;
end
Fnew = feval(FName,xnew,plist);
wnew = .5*(Fnew'*Fnew);
Jnew = feval(JName,xnew,plist);
function s = SineMercEarth(t) % s = SineMercEarth(t) % The sine of the Mercury-Sun-Earth angle at time t % Mercury location: xm = -11.9084 + 57.9117*cos(2*pi*t/87.97); ym = 56.6741*sin(2*pi*t/87.97); % Earth location: xe = -2.4987 + 149.6041*cos(2*pi*t/365.25); ye = 149.5832*sin(2*pi*t/365.25); s = (xm.*ye - xe.*ym)./(sqrt(xm.^2 + ym.^2).*sqrt(xe.^2 + ye.^2));
function s = DistMercEarth(t) % s = DistMercEarth(t) % The distance between Mercury and Earth at time t. % Mercury location: xm = -11.9084 + 57.9117*cos(2*pi*t/87.97); ym = 56.6741*sin(2*pi*t/87.97); % Earth location: xe = -2.4987 + 149.6041*cos(2*pi*t/365.25); ye = 149.5832*sin(2*pi*t/365.25); s = sqrt((xe-xm).^2 + (ye-ym).^2);
function A = TA(theta) % A = TA(theta) % % theta is a real number in the interval [-pi/2,pi/2]. % A is a negative multiple of the surface area of a tetrahedron % that is inscribed in the unit sphere. One vertex is at the north % pole (0,0,1) and the other three form an equilateral triangle in the plane % z = sin(theta). c = cos(theta); s = sin(theta); A = -c.*(sqrt((3*s-5).*(s-1)) + c);
function V = TV(theta) % V = TV(theta) % % theta real in interval [-pi/2,pi/2] % V is a negative multiple of the volume of a tetrahedron % that is inscribed in the unit sphere. One vertex is at the north % pole (0,0,1) and the other three form an equilateral triangle in the plane % z = sin(theta). s = sin(theta); V = (1-s.^2).*(s-1);
function E = TE(theta) % E = TE(theta) % % theta is a real number in the interval [-pi/2,pi/2] % E is a negative multiple of the total edge length of a tetrahedron % that is inscribed in the unit sphere. One vertex is at the north % pole (0,0,1) and the other three form an equilateral triangle in the plane % z = sin(theta). E = -(sqrt((1-sin(theta))) + sqrt(3/2)*cos(theta));
function E = TEDiff(theta) % E = TEDiff(theta) % % theta is a real number in the interval [-pi/2,pi/2] % % Let T be the tetrahedron with one vertex at the north pole (0.0,1) and the other % three forming an equilateral triangle in the plane z = sin(theta). E is difference % between a pole to base edge length and the length of a base edge. E = sqrt((1-sin(theta))) - sqrt(3/2)*cos(theta);
function pLoc = Orbit(t,planet,lineStyle)
% pLoc = Orbit(t,planet,lineStyle)
%
% t is a row vector and for k=1:length(t), pLoc.x(k) = x(t(k)) and
% pLoc.y(k) = y(t(k)) where
%
% x(tau) cos(phi) sin(phi) (planet.A-planet.P)/2 + ((planet.A+planet.P)/2)cos(tau)
% = *
% y(tau) -sin(phi) cos(phi) sqrt(planet.A*planet.P)sin(tau)
%
% If nargin==3 then the points are plotted with line style defined by the
% string lineStyle.
c = cos(t); s = sin(t);
x0 = ((planet.P-planet.A)/2) + ((planet.P+planet.A)/2)*c;
y0 = sqrt(planet.A*planet.P)*s;
cphi = cos(planet.phi); sphi = sin(planet.phi);
pLoc = struct('x',cphi*x0 + sphi*y0,'y',-sphi*x0 + cphi*y0);
if nargin==3
plot(pLoc.x,pLoc.y,lineStyle)
end
function d = Sep(t,planet1,planet2) % d = Sep(t,planet1,planet2) % % t is a 2-vector and planet1 and planet2 are orbit structures. % t(1) defines a planet1 orbit point, t(2) defines a planet2 orbit point, % and d is the distance between them. pLoc1 = Orbit(t(1),planet1); pLoc2 = Orbit(t(2),planet2); d = ((pLoc1.x-pLoc2.x)^2 + (pLoc1.y-pLoc2.y)^2)/2;
function g = gSep(t,planet1,planet2) % g = gSep(t,planet1,planet2) % The gradient of sep(t,planet1,planet2) with respect to 2-vector t. A1 = planet1.A; P1 = planet1.P; phi1 = planet1.phi; A2 = planet2.A; P2 = planet2.P; phi2 = planet2.phi; alfa1 = (P1-A1)/2; beta1 = (P1+A1)/2; gamma1 = sqrt(P1*A1); alfa2 = (P2-A2)/2; beta2 = (P2+A2)/2; gamma2 = sqrt(P2*A2); s1 = sin(t(1)); c1 = cos(t(1)); s2 = sin(t(2)); c2 = cos(t(2)); cphi1 = cos(phi1); sphi1 = sin(phi1); cphi2 = cos(phi2); sphi2 = sin(phi2); Rot1 = [cphi1 sphi1; -sphi1 cphi1]; Rot2 = [cphi2 sphi2; -sphi2 cphi2]; P1 = Rot1*[alfa1+beta1*c1;gamma1*s1]; P2 = Rot2*[alfa2+beta2*c2;gamma2*s2]; dP1 = Rot1*[-beta1*s1;gamma1*c1]; dP2 = Rot2*[-beta2*s2;gamma2*c2]; g = [-dP1';dP2']*(P2-P1);
function [g,H] = gHSep(t,planet1,planet2) % [g,H] = gHSep(t,planet1,planet2) % % t is a 2-vector and planet1 and planet2 are orbits. % g is the gradient of Sep(t,planet1,planet2) and H is the Hessian. A1 = planet1.A; P1 = planet1.P; phi1 = planet1.phi; A2 = planet2.A; P2 = planet2.P; phi2 = planet2.phi; alfa1 = (P1-A1)/2; beta1 = (P1+A1)/2; gamma1 = sqrt(P1*A1); alfa2 = (P2-A2)/2; beta2 = (P2+A2)/2; gamma2 = sqrt(P2*A2); s1 = sin(t(1)); c1 = cos(t(1)); s2 = sin(t(2)); c2 = cos(t(2)); cphi1 = cos(phi1); sphi1 = sin(phi1); cphi2 = cos(phi2); sphi2 = sin(phi2); Rot1 = [cphi1 sphi1; -sphi1 cphi1]; Rot2 = [cphi2 sphi2; -sphi2 cphi2]; P1 = Rot1*[alfa1+beta1*c1;gamma1*s1]; P2 = Rot2*[alfa2+beta2*c2;gamma2*s2]; dP1 = Rot1*[-beta1*s1;gamma1*c1]; dP2 = Rot2*[-beta2*s2;gamma2*c2]; g = [-dP1';dP2']*(P2-P1); ddP1 = Rot1*[-beta1*c1;-gamma1*s1]; ddP2 = Rot2*[-beta2*c2;-gamma2*s2]; H = zeros(2,2); H(1,1) = (P1(1)-P2(1))*ddP1(1) + (P1(2)-P2(2))*ddP1(2) + ... dP1(1)^2 + dP1(2)^2; H(2,2) = (P2(1)-P1(1))*ddP2(1) + (P2(2)-P1(2))*ddP2(2) + ... dP2(1)^2 + dP2(2)^2; H(1,2) = -dP1(1)*dP2(1) - dP1(2)*dP2(2); H(2,1) = H(1,2);
function z = SepV(t,planet1,planet2) % z = SepV(t,planet1,planet2) % The vector from the t(1) point on the planet1 orbit % to the t(2) point on the planet2 orbit. pt1 = Orbit(t(1),planet1); pt2 = Orbit(t(2),planet2); z = [pt2.x-pt1.x;pt2.y-pt1.y];
function J = JSepV(t,planet1,planet2) % J = JSepV(t,planet1,planet2) % J is the Jacobian of sepV(t,planet1,planet2). A1 = planet1.A; P1 = planet1.P; phi1 = planet1.phi; A2 = planet2.A; P2 = planet2.P; phi2 = planet2.phi; s1 = sin(t(1)); c1 = cos(t(1)); s2 = sin(t(2)); c2 = cos(t(2)); beta1 = (P1+A1)/2; gamma1 = sqrt(P1*A1); beta2 = (P2+A2)/2; gamma2 = sqrt(P2*A2); cphi1 = cos(phi1); sphi1 = sin(phi1); cphi2 = cos(phi2); sphi2 = sin(phi2); Rot1 = [cphi1 sphi1; - sphi1 cphi1]; Rot2 = [cphi2 sphi2; - sphi2 cphi2]; J = [ Rot1*[beta1*s1;-gamma1*c1] Rot2*[-beta2*s2;gamma2*c2]];
function f = rho(z,plist) % f = rho(z,plist) % z is a 2-vector whose components are the orbit parameters A and P % respectively. plist is a 2-column matrix. plist(i,1) is the i-th % observed radius vector length and plist(i,2) is the cosine of the % corresponding observed polar angle. A = z(1); P = z(2); r = plist(:,1); c = plist(:,2); f = r - (2*A*P)./(P*(1-c) + A*(1+c));
function J = Jrho(z,plist) % J = Jrho(z,plist) % J is the Jacobian of the function rho at z. A = z(1); P = z(2); c = plist(:,2); denom = (P*(1-c) + A*(1+c)).^2; J = -[ (2*P^2)*(1-c)./denom (2*A^2)*(1+c)./denom];