function [fp,f] = ksphere(Am, Ap, Fm0, Fp0, Iy , Ix, mu, looporder) %Updated 13 November 2014 %Computes the constant curvature surfaces in the 3-sphere using % the DPW (or generalized d'Alembert) method. The computation is %essentially the same as in the programme ksurf2, only a projection % to S^3 is taken, instead of applying the Sym-formula. The surface in %S^3 is the output f, and the output fp is a stereographic projection %to R^3, which is also plotted. %The loop group construction, with the same setup used here, is %described in arXiv:1301.5999 [math.DG] %"Constant Gaussian curvature surfaces in the 3-sphere via loop groups" %by David Brander, Jun-ichi Inoguchi, Shimpei Kobayashi. %All loops will be expanded up to order %'looporder' or the maximum order of the other input data in the loop parameter, %whichever is the greater. Increasing 'looporder' will result in greater %accuracy, but longer computation time. Also plots the surface using surf. %The loops described there %have to be untwisted before entering them into this program. This is %achieved by: twisted % [a(lambda), b(lambda); c(lambda), d(lambda)] % -> % [a(sqrt(lambda)), B_-(sqrt(lambda)); C_+(sqrt(lambda)), d(sqrt(lambda)] % untwisted % where B_-(lambda):= lambda^{-1} b(lambda) and %C_+(lambda):= lambda*c(lambda). %Author: David Brander; email: D.Brander@mat.dtu.dk. % Input: All loops are untwisted. % Ap = @(t) [A_-n(t),... A_n(t)] % where each A_i(t) is a 2 x 2 matrix valued function handle. %This represents the potential (A_-n lambda^{-1} + A_1 lambda + ... + A_n %lambda^n) dx. This corresponds to the "+" potential. % Similarly, Am = @(t) [B_-m (t) .... B_0 (t)], % represents the minus potential (B_-m lambda^(-m) + ... + B_0] dy. %"looporder" can be a positive integer "looporder=maxlooporder" %or a pair "looporder=[maxlooporder, minlooporder]" %(Default is minlooporder = floor(maxlooporder/2)). Loops %are approximated as polynomials of order minlooporder at the initial %point up to maxlooporder at the points furthest from the initial point. %Default for minlooporder is the 1 or the maximum order of Fp0 and Fm0. % All loops are converted into large square matrices, such that every pair of rows is a copy of the previous %pair, shifted to the right by 2 columns. Then matrix % multiplication corresponds exactly to loop multiplication, and the LU % decomposition corresponds to the birkhoff decomposition of loops. %Fp_0 and Fm_0 represent the initial conditions, and are Laurent %polynomial (untwisted) loops: Fp_0 = [F_-j , ...., F_j] % Fm_0 = [G_-k, ....., G_k], where F_i and G_i are constant 2x2 % matrices. Usually Fp_0 = Fm_0 = eye(2). %Ix = [x_init, x_step, leftpoints, rightpoints], where x_init is the initial x-value, %x_step is the step size between the x-values at which the %solution should be evaluated, and the integration is done in both directions, leftpoints to the left and rightpoints to the right of x_0. %If either of leftpoints or rightpoints is less than 3, then it will be %replaced by 0. %Iy is of a similar form. %fp, (the projection to R^3) is of the form [X; Y; Z], %where each of X,Y and Z are m x n matrices, %m is the number of x gridpoints 1+Ix(3)+Ix(4), and n the number of y %grid points 1+Iy(3)+Iy(4). %X is the matrix of x-values at the grid %points, Y the matrix of y values, and Z the matrix of z values. %This can be plotted with the command % surf(f(1:m, :), f(m+1: 2*m, :), f(2*m+1: 3*m, :)), % (that is, surf(X,Y,Z) ) which is also done automatically. %checkinput parses the input for consistency, converts the loops into %square matrix form and sets the minimum and maximum looporders. [xpoints,ypoints,bigAm, bigAp, bigFm0, bigFp0,minlooporder,maxlooporder]=checkinput(Am, Ap, Fm0, Fp0, Iy , Ix, looporder); ordermatrix = setlooporders(Ix,Iy,minlooporder,maxlooporder); %the orders for the LU decomposition and loop evaluation [Fplusinv,xerror]=invSintegrateA(bigAp, inv(bigFp0), Ix,2,1); disp(['x-axis: Left endpoint Error ', num2str(xerror(1),3),'. Right endpoint Error ', num2str(xerror(2),3),]); FpDims=size(Fplusinv); bsize=FpDims(1); [Fminus,yerror] = XSintegrateA(bigAm, bigFm0, Iy,2,1); disp(['y-axis: Left endpoint Error ', num2str(yerror(1),3),'. Right endpoint Error ', num2str(yerror(2),3),]); Y=zeros(2*xpoints,2*ypoints); ErrorMatrix=zeros(xpoints,ypoints); for j=1:xpoints factor1=looptranspose(Fplusinv(:,2*(j-1)+1:2*j),maxlooporder); for k=1:ypoints splitorder=ordermatrix(j,k); factor2=loop2SMat(Fminus(: ,(k-1)*bsize+2*maxlooporder+1:k*bsize-2*maxlooporder),splitorder,splitorder); Phi=factor1(:,2*(2*maxlooporder-splitorder)+1: 2*(2*maxlooporder+splitorder+1))*factor2; tosplit=loop2SMat(Phi,splitorder,splitorder); [~,right,~]=lu(tosplit,0.0000001); F = Fminus(: ,(k-1)*bsize+2*(2*maxlooporder-splitorder)+1:k*bsize-2*(2*maxlooporder-splitorder))/right; F1=rloopeval(F,1); Fmu=[1,mu;1/mu,1].*rloopeval(F,mu^2); f(2*j-1:2*j,2*k-1:2*k)=Fmu/F1; ErrorMatrix(j,k)=(norm(f(2*j-1:2*j,2*k-1:2*k)'*f(2*j-1:2*j,2*k-1:2*k)-eye(2)))/sqrt(2); end if mod(j,20)==1 disp(['Row ', num2str(j), ' Max Error ', num2str(max(ErrorMatrix(j:j,:)),2),'. Errors: ', num2str(ErrorMatrix(j:j,:),2)]); end end disp(' '); disp(['Max error:', num2str(max(max(ErrorMatrix)),2), '. Mean error: ', num2str(mean(mean(ErrorMatrix)),2), '. Corner Errors: ']); disp(num2str([ErrorMatrix(1,1), ErrorMatrix(1,ypoints);ErrorMatrix(xpoints,1), ErrorMatrix(xpoints,ypoints)],2)); fp = projplotdata(f); end %%%%%%%%%%%%%%%%%%%%%%%%%%%% sub functions %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Y = rloopeval(X,r) % Evaluates a loop at lambda = r. % Input: X = [x_-n .... x_0 x_1 .... x_n], where each x_i is 2x2 % r is a number. %Output: a sparse 2x2 matrix, the sum of x_i r^i n = (length(X)-2)/4; % the order of X in lambda Y= full(X(1:2, 2*n+1:2*n+2)); for count = 1 : n Y = Y + r^count*X(1:2, 2*n+1+2*count : 2*n+2+ 2*count) + r^(-count)*X(1:2, 2*n+1-2*count : 2*n + 2 - 2*count); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Y = projplotdata(f) % fheight=length(f(:,1:1)); xpoints=fheight/2; fwidth=length(f(1:1,:)); ypoints=fwidth/2; Y=zeros(3*xpoints, ypoints); for count=1:xpoints for j=1:ypoints Y(count,j) = -imag(f(2*count-1, 2*j-1))/(1+real(f(2*count-1, 2*j-1))); Y(xpoints + count, j) = -real(f(2*count-1, 2*j))/(1+real(f(2*count-1, 2*j-1))); Y(2*xpoints+count, j) = -imag(f(2*count-1, 2*j))/(1+real(f(2*count-1, 2*j-1))); end end surf(Y(1:xpoints, :), Y(xpoints+1: 2*xpoints, :), Y(2*xpoints+1: 3*xpoints, :),'EdgeAlpha',0.3,'FaceAlpha',1,'FaceLighting','gouraud','EdgeLighting','gouraud','EdgeColor','white','FaceColor',[0.05 0.4 0.8]) daspect([1 1 1]) grid off end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Finv,Error] = invSintegrateA(A, F0, Ix,m,lambda) %The integration along the y axis. %solves the equation dG/dt = -A(t) G %The values of Finv are returned as loops in column form. % Estimates the "error" at the endpoints by checking if the determinant is % 1. tempIy3=Ix(3); tempIy4=Ix(4); if Ix(3)<2, tempIy3=2; end if Ix(4)<2, tempIy4=2; end stepy= Ix(2); intervalyr=zeros(1,tempIy4+1); intervalyl=zeros(1,tempIy3+1); for j = 1 : tempIy4+1 intervalyr(j)= Ix(1) + (j-1)*stepy; end for j = 1 : tempIy3+1 intervalyl(j)= Ix(1) - (j-1)*stepy; end bsize=length(F0); n =(bsize/m-1)/4; % the order of A in lambda initF=reshape(F0(:, m*2*n+1:m*2*n+m),m*bsize,1); testfn = @(t,y) reshape((-A(t))*sparse(reshape(y,bsize,m)), m*bsize,1); if Ix(3)>0 [~,soll] = ode45(testfn ,intervalyl,initF); G1=sparse(reshape(flipud(soll(2:Ix(3)+1, :)).', bsize,m*(Ix(3)))); else G1=[]; end if Ix(4)>0 [~,solr] = ode45(testfn ,intervalyr,initF); G2=sparse(reshape(solr(2:Ix(4)+1, :).', bsize,m*(Ix(4)))); else G2=[]; end Finv=[G1, F0(:, m*2*n+1:m*2*n+m), G2]; Error=zeros(2,1); Error(1)=abs(1-det(rloopeval(looptranspose(Finv(:,1:m),n),lambda))); Error(2)=abs(1-det(rloopeval(looptranspose(Finv(:,m*(Ix(3)+Ix(4))+1:m*(Ix(3)+Ix(4)+1)),n),lambda))); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [F,Error] = XSintegrateA(A, F0, Ix,m,lambda) %solves the equation dF/dt = F A, with initial condition F(t0) = F0, with t %in the specified interval. % Estimates the "error" at the endpoints by checking if the determinant is % 1. tempIx3=Ix(3); tempIx4=Ix(4); if Ix(3)<2, tempIx3=2; end if Ix(4)<2, tempIx4=2; end stepy= Ix(2); intervalyl=zeros(1, tempIx3+1); intervalyr=zeros(1,tempIx4+1); for count = 1 : tempIx4+1 intervalyr(count)= Ix(1) + (count-1)*stepy; end for count=1: tempIx3+1 intervalyl(count)= Ix(1) - (count-1)*stepy; end bsize=length(F0); n =(bsize/m-1)/4; % the order of X in lambda initF=reshape(F0(m*2*n+1:m*2*n+m,:),m*bsize,1); y=zeros(m*bsize,1); testfn = @(t,y) reshape(reshape(y,m,bsize)*A(t), m*bsize,1); [~,sol1] = ode45(testfn ,intervalyl,initF); [~,sol2] = ode45(testfn ,intervalyr,initF); F=[reshape(flipud(sol1(2:Ix(3)+1, :)).', m,bsize*(Ix(3))), reshape(sol2(1:Ix(4)+1, :).', m,bsize*(Ix(4)+1))]; Error=zeros(2,1); Error(1)=abs(1-det(rloopeval(F(:,1:bsize),lambda))); Error(2)=abs(1-det(rloopeval(F(:,bsize*(Ix(3)+Ix(4))+1:bsize*(Ix(3)+Ix(4)+1)),lambda))); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function V=looptranspose(U,d) %input: U is a loop given as a block colum matrix % U= [U_(n); ... ; U_(-n)] % output: V is in row form V= [0,...,0,V_(-n+d), ... V_(n-d),0,...0] Udims=size(U); len=Udims(1); m=Udims(2); K=len/m; %the number of coefficients of U V=sparse(zeros(m,m*K)); for j=d:K-1-d V(:,m*j+1:m*j+m)=U(len-m*(j+1)+1:len-m*(j+1)+m,:); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Xmat = loop2matrix(X,n,k) %Converts a loop into a big square matrix %Input: a matrix X=[x_-nold .... x_0 x_1 .... x_nold] representing a Laurent %polynomial loop of order nold %The input loop does not need to be the same order as the output loop %x_-nold lambda^-nold + ... + x_0 + x_nold lambda^nold, where each x_i is an mxm matrix %Output: Xmat is a sparse m(4n+1) x m(4n+1) square matrix, which looks something %like [x_0 ... x_n 0 0 ... 0; % x_-1 x_0 ....x_n 0 ...; % x_-2 x_-1 x_0 .......; ..... etc % Multiplying Xmat with Ymat corresponds to multiplying the loops X and Y. Xdims=size(X); m=Xdims(1); nold =(Xdims(2)/m-1)/2; % the order of the input loop X in lambda k=min(nold,k); Row=X(1:m, m*(nold-k)+1: m*(nold+k+1)); Xmat=sparse(zeros(m*(4*n+1))); for j=0:n Xmat(m*n+m*j+1:m*n+m*j+m, m*(n-k+j)+1:m*(n-k+j+2*k+1))=Row; end Xmat(1:m*n, 1:m*(2*n+1))=Xmat(m*n+1:2*m*n,m*n+1:m*3*n+m); Xmat(m*(2*n+1)+1:m*(3*n+1), m*(n)+1:m*(4*n+1))=Xmat(m*n+m+1:m*(2*n+1),1:m*(3*n+1)); Xmat(m*(3*n+1)+1:m*(4*n+1), m*(2*n+1)+1:(4*n+1)*m)=Xmat(m*n+1:2*m*n,1:2*m*n); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Xmat = loop2SMat(X,n,k) %Converts a loop into a big square matrix %Input: a matrix X=[x_-nold .... x_0 x_1 .... x_nold] representing a Laurent %polynomial loop of order nold %The input loop does not need to be the same order as the output loop %x_-nold lambda^-nold + ... + x_0 + x_nold lambda^nold, where each x_i is an mxm matrix %Output: Xmat is a sparse m(4n+1) x m(4n+1) square matrix, which looks something %like [x_0 ... x_n 0 0 ... 0; % x_-1 x_0 ....x_n 0 ...; % x_-2 x_-1 x_0 .......; ..... etc % Multiplying Xmat with Ymat corresponds to multiplying the loops X and Y. Xdims=size(X); m=Xdims(1); nold =(Xdims(2)/m-1)/2; % the order of the input loop X in lambda k=min(nold,k); Xmat=zeros(m*(2*n+1)); Xmat(m*n+1:m*(n+1),m*(n-k)+1:m*(n+1+k))=full(X(1:m, m*(nold-k)+1: m*(nold+k+1))); for j=0:n-1 Xmat(m*j+1:m*j+m,1:m*(n+1+j))=Xmat(m*n+1:m*(n+1),m*(n-j)+1:m*(2*n+1)); Xmat(m*(2*n-j)+1:m*(2*n-j+1),m*(n-j)+1:m*(2*n+1))=Xmat(m*n+1:m*(n+1),1:m*(n+1+j)); end Xmat=sparse(Xmat); end %%%%%%%%%%%%%%%%%%%%%%%%%% function orders = setlooporders(Ix,Iy,minlooporder,maxlooporder) %creates a matrix of looporders to be used in the LU decomposition %to speed things up. orders=zeros(Ix(3)+Ix(4)+1,Iy(3)+Iy(4)+1); center=[Ix(3)+1,Iy(3)+1]; %the coordinates of the center of the matrix cornerdistances=[norm(center-[1,1]), norm(center-[1,Iy(3)+Iy(4)+1]), norm(center-[Ix(3)+Ix(4)+1,1]), norm(center-[Ix(3)+Ix(4)+1,Iy(3)+Iy(4)+1])]; for m=1:(Ix(3)+Ix(4)+1) for k=1:Iy(3)+Iy(4)+1 orders(m,k)=min(ceil(minlooporder+(maxlooporder-minlooporder)*norm(center-[m,k])/max(cornerdistances)),maxlooporder); end end end %%%%%%%%%%%%%%%%%%%%%%%%%% function [xpoints,ypoints,bigAm, bigAp, bigFm0, bigFp0,minlooporder,maxlooporder]=checkinput(Am, Ap, Fm0, Fp0, Iy , Ix, looporder) %Checks the input data for consistency with expected input. %Converts the intial data, in terms of simple untwisted loops, into %untwisted loops represented by big square matrices, expanded up to order %'looporder' or the maximum order of the other data in the loop parameter, %whichever is the greater. % Output: Each of bigAp and bigAm is a handle for a function of the form = @(t) M(t) % where M is a 2(4n+1) x 2(4n+1) square matrix valued function of t % which looks something like % [x_0 ... x_n 0 0 ... 0; % x_-1 x_0 ....x_n 0 ...; % x_-2 x_-1 x_0 .......; ..... etc, each x_i a 2x2 matrix % bigFp0 and bigFm0 are square matrices of initial values, of the same size %and similar structure to A. assert((Ix(3)>=0)&(Ix(4)>=0)&(Ix(3)+Ix(4)>0),'Ix(3) and Ix(4) must be non-negative, and at least one positive'); assert((Iy(3)>=0)&(Iy(4)>=0)&(Iy(3)+Iy(4)>0),'Iy(3) and Iy(4) must be non-negative, and at least one positive'); xpoints=1+Ix(3)+Ix(4); ypoints=1+Iy(3)+Iy(4); Fm0Dims=size(Fm0); Fp0Dims=size(Fp0); AmDims=size(Am(Iy(1))); ApDims=size(Ap(Ix(1))); m=Fm0Dims(1); assert((Fm0Dims(1)==2)&(Fp0Dims(1)==2)&(AmDims(1)==2)&(ApDims(1)==2),'All input loops must be matrices with two rows'); Fm0order = (Fm0Dims(2)-m)/(2*m); Fp0order = (Fp0Dims(2)-m)/(2*m); AmOrder = (AmDims(2)-m)/(2*m); ApOrder = (ApDims(2)-m)/(2*m); assert(norm(floor([Fm0order, Fp0order, AmOrder,ApOrder])-[Fm0order, Fp0order, AmOrder,ApOrder])==0,... 'Dimensions of input matrices and function handles are not consistent for loops'); loDims=size(looporder); assert((loDims(2)<=2)&(norm(floor(looporder)-looporder)==0)&(norm(looporder-abs(looporder))==0),... 'looporder must be a non-negative integer or a 1x2 non-negative integer matrix '); if loDims(2)==2 maxlooporder=max([looporder(1),Fm0order,Fp0order,AmOrder,ApOrder]); minlooporder=looporder(2); else maxlooporder=max([looporder,Fm0order,Fp0order,AmOrder,ApOrder]); minlooporder=max([floor(maxlooporder/2),Fm0order,Fp0order]); end %convert the loops into square matrix form: bigAm=@(t)loop2matrix(Am(t),maxlooporder,AmOrder); bigAp=@(t)loop2matrix(Ap(t),maxlooporder,ApOrder); bigFm0= loop2matrix(Fm0,maxlooporder,Fm0order); bigFp0= loop2matrix(Fp0,maxlooporder,Fp0order); end