Added chebfun library

Author: rahulvigneswaran

chebfun/+cheb/README.txt

@@ -0,0 +1,17 @@
+The purpose of this package is to provide extra capabilities not found in the
+core Chebfun classes and functions. To be added to this folder, a function must
+have the following properties:
+
+1. No other functions depend on it.
+
+2. It depends only on built-in MATLAB and core Chebfun.
+
+3. It is well documented, following MATLAB conventions for the help
+   documentation block. In particular, it should provide at least one example
+   of usage that will be printed out by MATLAB help.
+
+4. It has a comment with the author's name and email, and the month and year
+   of original authorship.
+
+5. It has a notice assigning copyright to The University of Oxford and The
+   Chebfun Developers.

chebfun/+cheb/bernoulli.m

@@ -0,0 +1,23 @@
+function B = bernoulli(N)
+%BERNOULLI   Bernoulli polynomials as chebfuns.
+%   B = BERNOULLI(N) returns a quasimatrix of the first N+1 Bernoulli
+%   polynomials on [0,1].
+%
+%   Example (Bernolli numbers):
+%      B = cheb.bernoulli(4);
+%      format rat, B(0,:)
+
+% Copyright 2017 by The University of Oxford and The Chebfun Developers.
+% See http://www.chebfun.org/ for Chebfun information.
+
+x = chebfun('x', [0,1]);
+
+% Initalize the constant 1 as the first Bernoulli polynomial.
+B = 0*x + 1;
+% Compute the requested polynomials.
+for j = 1:N
+    B(:,j+1) = j*cumsum(B(:,j));
+    B(:,j+1) = B(:,j+1) - sum(B(:,j+1));
+end
+
+end

chebfun/+cheb/bspline.m

@@ -0,0 +1,23 @@
+function B = bspline(m)
+%BSPLINE   B-splines as chebfuns.
+%   BSPLINE(M) returns the first M B-splines, starting from a constant value 1
+%   on [-1/2,1/2]. Each sucessive B-spline is smoother by one derivative and
+%   supported on an interval larger by 1/2 in each direction.
+%
+%   The returned result is a chebmatrix.
+%
+%   Example:
+%      B = cheb.bspline(4); hold on
+%      for k=1:4, plot(B{k}), end
+%      axis([-0.1 1.1 -2 2])
+
+% Copyright 2017 by The University of Oxford and The Chebfun Developers.
+% See http://www.chebfun.org/ for Chebfun information.
+
+s = chebfun(1, [-.5 .5]);
+B = chebmatrix(s);
+for k = 1:m-1
+    B(k+1) = conv(B{k}, s);
+end
+
+end

chebfun/+cheb/choreo.m

@@ -0,0 +1,187 @@
+function choreo
+%CHOREO   Compute planar choreographies of the n-body problem.
+%    CHEB.CHOREO computes a planar choreography using hand-drawn initial
+%    guesses. At the end of the computation, the user is asked to press <1> to
+%    display the motion of the planets. To stop the program, press <CTRL>-<C>.
+%
+% Planar choreographies are periodic solutions of the n-body problem in which
+% the bodies share a single orbit. This orbit can be fixed or rotating with some
+% angular velocity relative to an inertial reference frame.
+%
+% The algorithm uses trigonometric interpolation and quasi-Newton methods.
+% See [1] for details.
+%
+% Example: at the prompt, specify 5 bodies, angular rotation 0 or 1.2.
+% Then either draw a curve (if your machine has imfreehand) or click
+% in 10 or 15 points (if it doesn't) roughly along a figure-8.  Type
+% <Enter>, and in a few seconds you will see a choreograpy.
+%
+% [1] H. Montanelli and N. I. Gushterov, Computing planar and spherical
+% choreographies, SIAM Journal on Applied Dynamical Systems 15 (2016),
+
+% Copyright 2017 by The University of Oxford and The Chebfun Developers.
+% See http://www.chebfun.org/ for Chebfun information.
+
+close all
+LW = 'linewidth'; MS = 'markersize'; FS = 'fontsize';
+lw = 2; ms = 30; fs = 18;
+format long, format compact
+n = input('How many bodies? (e.g. 5) ');
+w = input('Angular velocity? (either 0 or a nonzero noninteger) ');
+k = 15;
+N = k*n;
+dom = [0 2*pi];
+
+while 1
+    
+    % Hand-drawn initial guess:
+    xmax = 2.5; ymax = 2.5;
+    clf, hold off, axis equal, axis([-xmax xmax -ymax ymax]), box on
+    set(gca,FS,fs), title('Draw a curve!')
+    if ( exist('imfreehand') == 2 )
+        % Try to use IMFREEHAND (need the image processing toolbox):
+        h = imfreehand();
+        z = getPosition(h);
+        delete(h)
+        z = z(:,1) + 1i*z(:,2);
+        q0 = chebfun(z,dom,'trig');
+        q0 = chebfun(q0,dom,N,'trig');
+        c0 = trigcoeffs(q0);
+        c0(1+floor(N/2)) = 0;
+        q0 = chebfun(c0,dom,'coeffs','trig');
+    else
+        % Otherwise, use GINPUT:
+        x = []; y = []; button = 1;
+        disp('Input points with mouse, press <enter> for final point.')
+        while ( button == 1 )
+            [xx,yy,button] = ginput(1);
+            x = [x; xx]; y = [y; yy]; %#ok<*AGROW>
+            hold on, plot(xx,yy,'xb',MS,ms/2), drawnow
+        end
+        z = x + 1i*y;
+        q0 = chebfun(z,dom,'trig');
+        q0 = chebfun(q0,dom,N,'trig');
+    end
+    
+    % Solve the problem:
+    hold on, plot(q0,'.-b',LW,lw), drawnow
+    c0 = trigcoeffs(q0);
+    c0 = [real(c0);imag(c0)];
+    [A0,G0] = actiongradeval(c0,n,w);
+    fprintf('\nInitial action: %.6f\n',A0)
+    options = optimoptions('fminunc');
+    options.Algorithm = 'quasi-newton';
+    options.HessUpdate = 'bfgs';
+    options.GradObj = 'on';
+    options.Display = 'off';
+    [c,A,~, ~,G] = fminunc(@(x)actiongradeval(x,n,w),c0,options);
+    fprintf('Action after optimization: %.6f\n',A)
+    fprintf('Norm of the gradient: %.3e\n',norm(G)/norm(G0))
+    
+    % Plot the result:
+    c = c(1:N) + 1i*c(N+1:2*N);
+    c(1+floor(N/2)) = 0;
+    q = chebfun(c,dom,'coeffs','trig');
+    hold on, plot(q,'r',LW,lw)
+    hold on, plot(q(2*pi*(0:n-1)/n),'.k',MS,ms), title('')
+    s = input('Want to see the planets dancing? (Yes=<1>, No=<0>) ');
+    if ( s == 1 )
+        hold off, plotonplane(q,n,w), pause
+    end
+    
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [A,G] = actiongradeval(c,n,w)
+%ACTIONGRADEVAL  Compute the action and its gradient in the plane.
+
+% Set up:
+if ( nargin < 3 )
+    w = 0;
+end
+N = length(c)/2;
+u = c(1:N);
+v = c(N+1:end);
+c = u + 1i*v;
+if ( mod(N, 2) == 0 )
+    k = (-N/2:N/2-1)';
+else
+    k = (-(N-1)/2:(N-1)/2)';
+end
+[t,s] = trigpts(N,[0 2*pi]);
+
+% Evaluate action A:
+c = bsxfun(@times,exp(1i*k*(0:n-1)*2*pi/n),c(:,1));
+vals = ifft(N*c);
+dc = 1i*k.*c(:,1);
+if ( mod(N, 2) == 0 )
+    dc(1,1) = 0;
+end
+dvals = ifft(N*dc);
+K = abs(dvals+1i*w*vals(:,1)).^2;
+U = 0;
+for i = 2:n
+    U = U - 1./abs(vals(:,1)-vals(:,i));
+end
+A = n/2*s*(K-U);
+
+if ( nargout > 1 )
+    
+    % Initialize gradient G:
+    G = zeros(2*N,1);
+    
+    % Loop over the bodies for the AU contribution:
+    for j = 1:n-1
+        a = bsxfun(@times,1-cos(k'*j*2*pi/n),cos(t*k')) + ...
+            bsxfun(@times,sin(k'*j*2*pi/n),sin(t*k'));
+        b = bsxfun(@times,-1+cos(k'*j*2*pi/n),sin(t*k')) + ...
+            bsxfun(@times,sin(k'*j*2*pi/n),cos(t*k'));
+        f = (a*u + b*v).^2 + (-b*u + a*v).^2;
+        df = 2*(bsxfun(@times,a*u + b*v,a) + bsxfun(@times,-b*u + a*v,-b));
+        G(1:N) = G(1:N) - n/4*bsxfun(@rdivide,df,f.^(3/2))'*s';
+        df = 2*(bsxfun(@times,a*u + b*v, b) + bsxfun(@times,-b*u + a*v,a));
+        G(N+1:2*N) = G(N+1:2*N) - n/4*bsxfun(@rdivide,df,f.^(3/2))'*s';
+    end
+    
+    % Add the AK contribution:
+    if ( mod(N, 2) == 0 )
+        G = G + 2*pi*n*[w;(k(2:end)+w);w;(k(2:end)+w)].^2.*[u;v];
+    else
+        G = G + 2*pi*n*[(k+w);(k+w)].^2.*[u;v];
+    end
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function plotonplane(q,n,w)
+%PLOTONPLANE   Plot a planar choreography.
+
+vscale = max(max(real(q)),max(imag(q)));
+xmax = 1.3*vscale; ymax = xmax;
+LW = 'linewidth'; MS = 'markersize'; FS = 'fontsize';
+lw = 2; ms = 30; fs = 18;
+dom = [0 2*pi]; T = 2*pi; dt = .1;
+xStars = 2*xmax*rand(250,1)-xmax;
+yStars = 2*ymax*rand(250,1)-ymax;
+for t = dt:dt:10*T
+    clf, hold off
+    fill(xmax*[-1 1 1 -1 -1],ymax*[-1 -1 1 1 -1],'k')
+    axis equal, axis([-xmax xmax -ymax ymax]), box on
+    set(gca,FS,fs)
+    hold on, plot(xStars,yStars,'.w',MS,4)
+    if ( w ~= 0 )
+        q = chebfun(@(t)exp(1i*w*dt).*q(t),dom,length(q),'trig');
+    end
+    if ( t < T )
+        hold on, plot(q,'r',LW,lw)
+    end
+    for j = 1:n
+        hold on, plot(q(t+2*pi*j/n),'.',MS,ms)
+    end
+    drawnow
+end
+
+end