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];