No BSD License  

Highlights from
Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

image thumbnail

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

by

Howard Wilson

 

14 Oct 2002 (Updated )

Companion Software (amamhlib)

brachist
function brachist      
% Example: brachist
% ~~~~~~~~~~~~~~~~~
% This program determines the shape of a 
% smooth curve down which a particle can slide 
% in minimum possible time. The analysis 
% employs a piecewise cubic spline to define
% the curve by interpolation using a fixed set 
% of base point positions. The curve shape 
% becomes a function of the heights chosen at 
% the base points. These heights are determined 
% by computing the descent time as a function
% of the heights and minimizing the descent 
% time by use of an optimization program. The 
% Nelder and Mead unconstrained search 
% procedure is used to compute the minimum.
%
% User m functions called:  
%      chbpts, brfaltim, fltim, gcquad, 
%      bracifun, splined

global cbp cwf cofs n xc yc a b b_over_a ...
       grav nparts nquad nfcls 

fprintf(...
'\nBRACHISTOCHRONE DETERMINATION BY NONLINEAR');
fprintf('\n         OPTIMIZATION SEARCH \n');
fprintf(['\nPlease wait. The ',...
  'calculation takes a while.\n']);

% Initialize
a=30; b=10; grav=32.2; nparts=1; nquad=50; 
tol=1e-4; n=6; b_over_a = b/a;

[dummy,cbp,cwf]=gcquad('',0,1,nquad,nparts);
xc=chbpts(0,1,n); xc=xc(:);
y0=5*sin(pi*xc); xc=[0;xc;1];

% Calculate results from the exact solution
[texact,xexact,yexact]=brfaltim(a,b,grav,100);

% Perform minimization search for 
% approximate solution
opts=optimset('tolx',tol,'tolfun',tol);
[yfmin,fmin,flag,outp] =...
                   fminsearch(@fltim,y0,opts);

% Evaluate final position and approximate 
% descent time
Xfmin=xc; Yfmin=Xfmin+[0;yfmin(:);0];
% tfmin=a/sqrt(2*grav*b)*fltim(yfmin(:));
tfmin=a/sqrt(2*grav*b)*fmin;
nfcls=1+outp.funcCount;

% Summary of calculations
fprintf('\nBrachistochrone Summary');
fprintf('\n-----------------------');
fprintf('\nNumber of function calls:   ');
fprintf('%g',nfcls);
fprintf('\nDescent time:   ');
fprintf('%g',tfmin), fprintf('\n')

% Plot results comparing the approximate 
% and exact solutions
xplot=linspace(0,1,100); close
yplot=spline(Xfmin,Yfmin,xplot); 
plot(xplot,-yplot,'-',Xfmin,-Yfmin,'o', ...
     xexact/a,-yexact/b,'--');
xlabel('x/a'); ylabel('y/b'); % grid
title(['Brachistochrone Curve for ', ...
       'a/b = ',num2str(a/b)]);
text(.5,-.1,   'Descent time    (secs)')
text(.5,-.175,['Approximate: ',num2str(tfmin)])
text(.5,-.25, ['Exact:       ',num2str(texact)]);
text(.5,-.325, ...
  sprintf('Error tolerance:    %g',tol));
legend('Approximate Curve', ...
       'Computed Points','Exact Curve',3); 
figure(gcf); 
print -deps brachist

%=============================================

function [tfall,xbrac,ybrac]=brfaltim ...
                             (a,b,grav,npts)
%
% 
% [tfall,xbrac,ybrac]=brfaltim(a,b,grav,npts)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 
% This function determines the descent time 
% and a set of points on the brachistochrone 
% curve passing through (0,0) and (a,b). 
% The curve is a cycloid expressible in 
% parametric form as
%
%    x=k*(th-sin(th)), 
%    y=k*(1-cos(th))    for 0<=th<=thf
%
% where thf is found by solving the equation
%
%    b/a=(1-cos(thf))/(thf-sin(thf)). 
%
% Once thf is known then k is found from 
%
%    k=a/(th-sin(th)). 
%
% The exact value of the descent time is given 
% by
%
%    tfall=sqrt(k/g)*thf
%  
% a,b   - final values of (x,y) on the curve
% grav  - the gravity constant
% npts  - the number of points computed on 
%         the curve
%
% tfall - the time required for a smooth 
%         particle to slide along the curve 
%         from (0,0) to (a,b)
% xbrac - x points on the curve with x 
%         increasing to the right
% ybrac - y points on the curve with y 
%         increasing downward           
%
% User m functions called: none
%----------------------------------------------

brfn=inline('cos(th)-1+cof*(th-sin(th))','th','cof');
 
ba=b/a; [th,fval,flag]=fzero(...
                 brfn,[.01,10],optimset('fzero'),ba);

k=a/(th-sin(th)); tfall=sqrt(k/grav)*th;
if nargin==4
  thvec=(0:npts-1)'*(th/(npts-1));
  xbrac=k*(thvec-sin(thvec)); 
  ybrac=k*(1-cos(thvec));
end

%=============================================

function x=chbpts(xmin,xmax,n)
%
% x=chbpts(xmin,xmax,n)
% ~~~~~~~~~~~~~~~~~~~~~
% Determine n points with Chebyshev spacing 
% between xmin and xmax.
%
% User m functions called:  none
%----------------------------------------------

x=(xmin+xmax)/2+((xmin-xmax)/2)* ...
  cos(pi/n*((0:n-1)'+.5));

%=============================================

function t=fltim(y)
%
% t=fltim(y)
% ~~~~~~~~~~
%
% This function evaluates the time descent 
% integral for a spline curve having heights 
% stored in y.
%
% y - vector defining the curve heights at 
%     interior points corresponding to base 
%     positions in xc
%
% t - the numerically integrated time descent 
%     integral evaluated by use of base points 
%     cbp and weight factors cwf passed as 
%     global variables
%
% User m functions called: splined
%----------------------------------------------

global xc cofs nparts bp wf nfcls cbp cwf ...
       b_over_a

nfcls=nfcls+1; x=cbp;

% Generate coefficients used in spline 
% interpolation
yc=[0;y(:);0];
y=spline(xc,yc,x); yp=splined(xc,yc,x);

% Evaluate the integrand 
f=(1+(b_over_a*(1+yp)).^2)./(x+y); f=sqrt(f);

% Evaluate the integral
t=cwf(:)'*f(:);

%==============================================

% function [val,bp,wf]=gcquad(func,xlow,...
%                    xhigh,nquad,mparts,varargin)
% See Appendix B

%==============================================

% function val=splined(xd,yd,x,if2)
% See Appendix B 

Contact us