diff --git a/documentation/continuationInstructions.mlx b/documentation/continuationInstructions.mlx index 4d7bf971..c8a82640 100644 Binary files a/documentation/continuationInstructions.mlx and b/documentation/continuationInstructions.mlx differ diff --git a/resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolveDAEDopri8.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/highFreqFilter.m.type.File.xml similarity index 100% rename from resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolveDAEDopri8.m.type.File.xml rename to resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/highFreqFilter.m.type.File.xml diff --git a/resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolve15s.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolve15s.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolve15s.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolveARK.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolveARK.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolveARK.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/HHBDetector.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/HHBDetector.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/HHBDetector.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/HHBLocator.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/HHBLocator.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/HHBLocator.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File.xml new file mode 100644 index 00000000..1c0844ef --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File.xml @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/1.type.DIR_SIGNIFIER.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/1.type.DIR_SIGNIFIER.xml new file mode 100644 index 00000000..1c0844ef --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/1.type.DIR_SIGNIFIER.xml @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/HHBDetector.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/HHBDetector.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/HHBDetector.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/HHBLocator.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/HHBLocator.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/HHBLocator.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/adjointContinuationMethodFast.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/adjointContinuationMethodFast.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/adjointContinuationMethodFast.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/adjointContinuationMethodShift.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/adjointContinuationMethodShift.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/adjointContinuationMethodShift.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/eigsJL.m.type.File.xml b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/eigsJL.m.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/eigsJL.m.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File.xml b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File.xml new file mode 100644 index 00000000..1c0844ef --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File.xml @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/1.type.DIR_SIGNIFIER.xml b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/1.type.DIR_SIGNIFIER.xml new file mode 100644 index 00000000..1c0844ef --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/1.type.DIR_SIGNIFIER.xml @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/dumbbellTimeEvolutionDemo.mlx.type.File.xml b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/dumbbellTimeEvolutionDemo.mlx.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/dumbbellTimeEvolutionDemo.mlx.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/heatOnDumbbell.mlx.type.File.xml b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/heatOnDumbbell.mlx.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/heatOnDumbbell.mlx.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/lineTimeEvolutionDemo.mlx.type.File.xml b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/lineTimeEvolutionDemo.mlx.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/lineTimeEvolutionDemo.mlx.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/lollipopTimeEvolutionDemo.mlx.type.File.xml b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/lollipopTimeEvolutionDemo.mlx.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/lollipopTimeEvolutionDemo.mlx.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/starTimeEvolutionDemo.mlx.type.File.xml b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/starTimeEvolutionDemo.mlx.type.File.xml new file mode 100644 index 00000000..80b5b161 --- /dev/null +++ b/resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File/starTimeEvolutionDemo.mlx.type.File.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/resources/project/Root.type.ProjectPath/d4264395-e014-41c5-8cf6-74b2dd45f9a9.type.Reference.xml b/resources/project/Root.type.ProjectPath/d4264395-e014-41c5-8cf6-74b2dd45f9a9.type.Reference.xml new file mode 100644 index 00000000..2bd95e8d --- /dev/null +++ b/resources/project/Root.type.ProjectPath/d4264395-e014-41c5-8cf6-74b2dd45f9a9.type.Reference.xml @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/resources/project/Root.type.ProjectPath/d6c4416b-d461-4901-ae3f-0511106fdf87.type.Reference.xml b/resources/project/Root.type.ProjectPath/d6c4416b-d461-4901-ae3f-0511106fdf87.type.Reference.xml deleted file mode 100644 index 623c4230..00000000 --- a/resources/project/Root.type.ProjectPath/d6c4416b-d461-4901-ae3f-0511106fdf87.type.Reference.xml +++ /dev/null @@ -1,2 +0,0 @@ - - \ No newline at end of file diff --git a/resources/project/Root.type.ProjectPath/e6a4ab36-b927-44c8-9c3b-d8c041cd946a.type.Reference.xml b/resources/project/Root.type.ProjectPath/e6a4ab36-b927-44c8-9c3b-d8c041cd946a.type.Reference.xml new file mode 100644 index 00000000..0c8be5c7 --- /dev/null +++ b/resources/project/Root.type.ProjectPath/e6a4ab36-b927-44c8-9c3b-d8c041cd946a.type.Reference.xml @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/source/@quantumGraph/addChebyshevCoordinates.m b/source/@quantumGraph/addChebyshevCoordinates.m index 1a727dbd..fc61e173 100644 --- a/source/@quantumGraph/addChebyshevCoordinates.m +++ b/source/@quantumGraph/addChebyshevCoordinates.m @@ -20,7 +20,7 @@ function addChebyshevCoordinates(G,nxVec,ny) for k=1:nEdges nx=nxVec(k); - G.qg.Edges.x{k} = chebptsSecondKind(nx)'; + G.qg.Edges.x{k} = LVec(k)*chebptsSecondKind(nx)'; G.qg.Edges.y{k} = nan(nxVec(k),ny); end G.qg.Nodes.y=nan(nNodes,1); \ No newline at end of file diff --git a/source/@quantumGraph/column2graph.m b/source/@quantumGraph/column2graph.m index 200ad067..c349ed92 100644 --- a/source/@quantumGraph/column2graph.m +++ b/source/@quantumGraph/column2graph.m @@ -3,7 +3,7 @@ function column2graph(G,col) % a quantum graph with the same structure as template [~,nxC,nxTot] = G.nx; -assert(length(col) == nxTot,'Column and template not compatible.') +% assert(length(col) == nxTot,'Column and template not compatible.') for k=1:numedges(G) G.qg.Edges.y{k}=col(nxC(k)+1:nxC(k+1)); diff --git a/source/@quantumGraph/constructMatricesChebyshev.m b/source/@quantumGraph/constructMatricesChebyshev.m index df133b56..1246d7c0 100644 --- a/source/@quantumGraph/constructMatricesChebyshev.m +++ b/source/@quantumGraph/constructMatricesChebyshev.m @@ -10,7 +10,7 @@ function constructMatricesChebyshev(G) nNodes = G.numnodes; % Number of nodes [n,nxC,nxTot] = G.nx; % A useful vector giving positions of final disc point of each edge L = G.L; % Vector of edge lengths -D1matrix = zeros(nxTot,nxTot,nEdges); +D1matrix = zeros(nxTot,nxTot); D2matrix = zeros(nxTot-2*nEdges, nxTot); % D2matrix will contain the blocked D2 matricies and no BCs B = zeros(nxTot,nxTot); % Initializes projection matrix @@ -31,7 +31,7 @@ function constructMatricesChebyshev(G) x = x + n(i-1); % Positions us so that we move past the previously built portion of the D2 part of the matrix end - D1matrix( (nxC(i)+1):nxC(i+1) , (nxC(i)+1):nxC(i+1) , i) = -D1; % Creates square D1 matrix for edge_i leaving zeros elsewhere + D1matrix( (nxC(i)+1):nxC(i+1) , (nxC(i)+1):nxC(i+1)) = D1; % Creates square D1 matrix D2matrix( (x+1):(x+M) , (nxC(i)+1):nxC(i+1) ) = D2; B( (x+1):(x+M) , (nxC(i)+1):nxC(i+1) ) = Project; @@ -49,19 +49,19 @@ function constructMatricesChebyshev(G) % Create the Vertex Condition Assignment Matrix VCAMat which maps the nonhomogeneous data % defined at the vertices to the correct row -VCAMat = zeros(nxTot,nNodes,nNodes); +VCAMat = zeros(nxTot,nNodes); BC = zeros(2*nEdges,nxTot); % Initializes space for Boundary Conditions row = 0; for j=1:nNodes % Loop over the nodes [fullDegree,~,~] = G.fullDegreeEtc(j); for k=1:fullDegree % Loop over the edges connected to the node - row = row+1; - + row = row + 1; + if k == 1 % At first entry of block, enforce either Dirichlet or flux condition & put a one in the right spot of VCAMat BC(row,:) = robinCondition(G,j,D1matrix); - VCAMat(row,j) = 1; - else % At remaining entries of block, enforce continuity condition + VCAMat(nxTot-2*nEdges+row,j) = 1; + else % At remaining entries of block enforce continuity condition e1 = G.ek(j,1); e2 = G.ek(j,k); BC(row,:) = e1 - e2; @@ -71,10 +71,13 @@ function constructMatricesChebyshev(G) end LMat = [D2matrix; BC]; +C = B; +C(nxTot-2*nEdges+1:end,:) = BC; - +G.derivative = -D1matrix; G.laplacianMatrix = LMat; G.weightMatrix = B; +G.weightMatrixWithBCs = C; G.vertexConditionAssignmentMatrix=VCAMat; diff --git a/source/@quantumGraph/constructMatricesUniform.m b/source/@quantumGraph/constructMatricesUniform.m index 1fa107a9..6df88c25 100644 --- a/source/@quantumGraph/constructMatricesUniform.m +++ b/source/@quantumGraph/constructMatricesUniform.m @@ -69,7 +69,7 @@ function constructMatricesUniform(G) [ghostPt,interiorPt] = G.getEndPoints(direction,edge); row=row+1; - if k ==1 % At first entry of block, enforce either Dirichlet or flux condition & put a one in the right spot of VCAMat + if k == 1 % At first entry of block, enforce either Dirichlet or flux condition & put a one in the right spot of VCAMat if G.isDirichlet(j) % Define Dirichlet Boundary Condition LMat(row,ghostPt)=1/2; LMat(row,interiorPt)=1/2; diff --git a/source/@quantumGraph/dotChebyshev.m b/source/@quantumGraph/dotChebyshev.m index 5029e415..8c549cac 100644 --- a/source/@quantumGraph/dotChebyshev.m +++ b/source/@quantumGraph/dotChebyshev.m @@ -1,4 +1,4 @@ -function z = qgdotCheb(Phi,u1col,u2col) +function z = dotChebyshev(Phi,u1col,u2col) % Note this does absolutely no checking to determine if two columns are % compatible with the template Phi. It just surges ahead and returns an % error if it doesn't work @@ -18,9 +18,14 @@ function z = edgeDotCheb(Phi,v1,v2,k) - - l = Phi.Edges.L(k); + f = v1.*v2; - z = clencurt(f,l); % Uses Clenshaw-Curtis quadrature to integrate + l = Phi.Edges.L(k); + w = clencurtWeights(length(f)); % Weights for Clenshaw-Curtis quadrature + z = l/2 * w' * f; % Clenshaw-Curtis quadrature + +% l = Phi.Edges.L(k); +% f = v1.*v2; +% z = clencurt(f,l); % Uses Clenshaw-Curtis quadrature to integrate end \ No newline at end of file diff --git a/source/@quantumGraph/eigs.m b/source/@quantumGraph/eigs.m index 7ccde3e3..0cacc4d2 100644 --- a/source/@quantumGraph/eigs.m +++ b/source/@quantumGraph/eigs.m @@ -9,13 +9,13 @@ vec = zeros(length(A),n); -if isempty(null(full(A))) % If A is well conditioned - [V,d] = eigs(A,B,n,'SM'); % Use regular eigs solver as it is good enough +if isempty(null(full(A))) % If A is well conditioned + [V,d] = eigs(A,B,n,'SM'); % Use regular eigs solver as it is good enough val = diag(real(d)); -else % If A is ill conditioned - Ashift = A - B; % Shift A using B so Ashift is not be singular - [V,d] = eigs(Ashift,B,n,'SM'); % Solves (-A + 1B)u = lambda Bu - val = diag(real(d)) + 1 ; % Shift evals back and neglect small imaginary component +else % If A is ill conditioned + Ashift = A - B; % Shift A using B so Ashift is not be singular + [V,d] = eigs(Ashift,B,n,'SM'); % Solves (-A + 1B)u = lambda Bu + val = diag(real(d)) + 1 ; % Shift evals back and neglect small imaginary component end k=1; diff --git a/source/@quantumGraph/ek.m b/source/@quantumGraph/ek.m index 85ac5f3b..6a42996e 100644 --- a/source/@quantumGraph/ek.m +++ b/source/@quantumGraph/ek.m @@ -1,9 +1,12 @@ function ek = ek(G,node,k) +% Given a node (node) and kth edge affiliated with node, ek returns the standard +% basis vector coresponding to the discretization point of edge_k based on +% edge's orientation with respect to node. [~,nxC,nxTot] = G.nx; I = eye(nxTot); - [fullDegree,inOrOut,adjacentEdges] = G.fullDegreeEtc(node); + [~,inOrOut,adjacentEdges] = G.fullDegreeEtc(node); - if inOrOut(k) == 1 % If edge is incoming... + if inOrOut(k) == 1 % If edge is incoming... ek = I(nxC(adjacentEdges(k)+1),:); % ek gives us the last disc point on the edge else % If edges is outgoing... ek = I(nxC(adjacentEdges(k))+1,:); % ek gives us the first disc point on the edge diff --git a/source/@quantumGraph/isDirichlet.m b/source/@quantumGraph/isDirichlet.m index a2bb3416..f03961bf 100644 --- a/source/@quantumGraph/isDirichlet.m +++ b/source/@quantumGraph/isDirichlet.m @@ -3,5 +3,5 @@ if nargin==1 z = isnan(G.robinCoeff); else - z=isnan(G.robinCoeff(varargin{1})); + z = isnan(G.robinCoeff(varargin{1})); end \ No newline at end of file diff --git a/source/@quantumGraph/quantumGraph.m b/source/@quantumGraph/quantumGraph.m index 13d3cb55..61532a36 100644 --- a/source/@quantumGraph/quantumGraph.m +++ b/source/@quantumGraph/quantumGraph.m @@ -6,8 +6,10 @@ properties qg % the main quantum graph, a digraph discretization; + derivative; laplacianMatrix; weightMatrix; + weightMatrixWithBCs; vertexConditionAssignmentMatrix; explicitLaplacian; end diff --git a/source/@quantumGraph/robinCondition.m b/source/@quantumGraph/robinCondition.m index dec8dc99..47f897c4 100644 --- a/source/@quantumGraph/robinCondition.m +++ b/source/@quantumGraph/robinCondition.m @@ -1,20 +1,20 @@ function robin = robinCondition(G,node,D1matrix) -% Finds the Robin Condition for the given node +% Finds the discretized Robin Condition for the given node. alpha = G.robinCoeff; [fullDegree,inOrOut,adjacentEdges] = G.fullDegreeEtc(node); if G.isDirichlet(node) - robin = G.ek(node,adjacentEdges(1)); + robin = G.ek(node,1); else for k = 1:fullDegree c = inOrOut(k); w = G.Edges.Weight(adjacentEdges(k)); ek = G.ek(node,k); if k == 1 - robin = c*w*ek*D1matrix(:,:,adjacentEdges(1)) + c*alpha(node)*ek; + robin = c*w*ek*D1matrix + alpha(node)*ek; else - robin = robin + c*w*ek*D1matrix(:,:,adjacentEdges(k)); + robin = robin + c*w*ek*D1matrix; end end end diff --git a/source/@quantumGraph/secularDet.m b/source/@quantumGraph/secularDet.m index 68a6dea2..1fc4f121 100644 --- a/source/@quantumGraph/secularDet.m +++ b/source/@quantumGraph/secularDet.m @@ -28,7 +28,7 @@ for jprime =1:2*n % Builds S if dv==1 && jprime==jbar if j<=n && ~isnan(G.robinCoeff(G.EndNodes(j,2))) % Robin Boundary Conditions - alpha = G.robinCoeff(G.EndNodes(j,2)); + alpha = -G.robinCoeff(G.EndNodes(j,2)); S(jprime,j) = (1i*x+alpha)/(1i*x-alpha); if alpha == 0 ThetaR = ThetaR+pi; @@ -36,7 +36,7 @@ ThetaR = ThetaR + atan(alpha/x); end elseif jprime<=n && ~isnan(G.robinCoeff(G.EndNodes(jprime,1))) - alpha = G.robinCoeff(G.EndNodes(jprime,1)); + alpha = -G.robinCoeff(G.EndNodes(jprime,1)); S(jprime,j) = (1i*x+alpha)/(1i*x-alpha); if alpha == 0 ThetaR = ThetaR+pi; @@ -52,7 +52,7 @@ node = G.sharedNode(j,jprime); % Kirchhoff condition alpha = -G.robinCoeff(node); S(jprime,j) = 2/(dv-alpha/(1i*x))-1; - ThetaK = ThetaK + atan(alpha/(dv*x))/dv; % Divide by dv because the code loops through edges so we'll do this for each edge + ThetaK = ThetaK + atan(alpha/(dv*x))/dv; % Divide by dv because the code loops through edges so this is done for each edge elseif G.follows(j,jprime) node = G.sharedNode(j,jprime); % Kirchhoff condition alpha = -G.robinCoeff(node); diff --git a/source/@quantumGraph/timeEvolve15s.asv b/source/@quantumGraph/timeEvolve15s.asv new file mode 100644 index 00000000..bae396be --- /dev/null +++ b/source/@quantumGraph/timeEvolve15s.asv @@ -0,0 +1,30 @@ +function [t,u] = timeEvolve15s(G,f,u0,tend) +% Given a function f and initial condition u0, this function evolves the +% solution to f(t,x) = u_xx in time using a ode15s until t=tend using the +% given step size h. The result will be a matrix the ith row of u is the +% solution evaluated at ith value of tvec. + + if isUniform(G) + tol = 10-07; + else % Ensures the tolerence is set as small as possible + nx = G.Edges.nx(1); % given how accurately the elliptic problem can + if nx<=16 % be solved for a given number of discretization pts + tol = 10^-13; + elseif nx<=32 + tol = 10^-12; + else + tol = 10^-08; + end + end + + if lenth() + tspan = [0,tend]; + end + B = G.weightMatrix; + odeopt = odeset('mass', B, 'masssing', 'yes', 'AbsTol', tol, 'RelTol', tol); + + [t, u] = ode15s(f, tspan, u0, odeopt); + + u = transpose(u); + +end \ No newline at end of file diff --git a/source/@quantumGraph/timeEvolve15s.m b/source/@quantumGraph/timeEvolve15s.m new file mode 100644 index 00000000..cb17efb7 --- /dev/null +++ b/source/@quantumGraph/timeEvolve15s.m @@ -0,0 +1,32 @@ +function [t,u] = timeEvolve15s(G,f,u0,t) +% Given a function f and initial condition u0, this function evolves the +% solution to f(t,x) = u_xx in time using a ode15s until t=tend using the +% given step size h. The result will be a matrix the ith row of u is the +% solution evaluated at ith value of tvec. + + if isUniform(G) + tol = 10-07; + else % Ensures the tolerence is set as small as possible + nx = G.Edges.nx(1); % given how accurately the elliptic problem can + if nx<=16 % be solved for a given number of discretization pts + tol = 10^-13; + elseif nx<=32 + tol = 10^-12; + else + tol = 10^-08; + end + end + + if length(t) == 1 + tspan = [0,t]; + else + tspan = t; + end + B = G.weightMatrix; + odeopt = odeset('mass', B, 'masssing', 'yes', 'AbsTol', tol, 'RelTol', tol); + + [t, u] = ode15s(f, tspan, u0, odeopt); + + u = transpose(u); + +end \ No newline at end of file diff --git a/source/@quantumGraph/timeEvolveARK.asv b/source/@quantumGraph/timeEvolveARK.asv new file mode 100644 index 00000000..f728ce89 --- /dev/null +++ b/source/@quantumGraph/timeEvolveARK.asv @@ -0,0 +1,108 @@ +function [t,u] = timeEvolveARK(G,f,u0,h,tend) +% Given a graph G, function f and initial condition u0, this function +% evolves the solution to f(t,x) = u_xx in time using a 2nd, 4th or 8th +% order scheme until t=tend using the given step size h. The result will be +% a matrix the ith row of u is the solution evaluated at ith value of tvec. + +told = 0; % t_n +tnew = told; % t_n+1 +s = 3; % Number of b_i coeff's +j = 1; % Index for t and u +m = round(tend/h); % Number of t points +n = length(u0); % Number of x points + +nEdges = length(G.L); % Number of edges +nBCs = 2*nEdges; % Number of BC conditions +A = full(G.laplacianMatrix); % The second spatial derivative operator +B = full(G.weightMatrix); % The rectangular collocation or weight matrix +C = [B(1:(n-nBCs),:); A((n-nBCs+1):n,:)]; % Weight matrix diagonal with boundary data in bottom rows + + + +if size(u0,1)==1 % Corrects orientation of u0 + u0=u0'; +end + +u = zeros(n,m); +t = zeros(1,m); +u(:,1) = u0; +uold = u(:,1); % u_n + +x = 1.06858; +a = [x 0 0; 1/2-x x 0; 2*x 1-4*x x]; +b1 = 1/(6*(1-2*x)^2); +b2 = (3*(1-2*x)^2 - 1) / (3*(1-2*x)^2); +b3 = b1; +b = [b1 b2 b3]; +c = zeros(s,1); + +for i=1:s + c(i) = sum(a(i,:)); +end + + +while tnew < tend + told = tnew; + + k1guess = f(told , uold); + k1 = fsolve(@(x) x - f(told + c(1)*h, uold + a(1,1)*h*(C\x)), k1guess); + k2 = fsolve(@(x) x - f(told + c(2)*h, uold + a(2,1)*h*(C\k1) + a(2,2)*h*(C\x)), k1); + k3 = fsolve(@(x) x - f(told + c(3)*h, uold + a(3,1)*h*(C\k1) + a(3,2)*h*(C\k2) + a(3,3)*h*(C\x)), k2); + + k = [transpose(k1); transpose(k2); transpose(k3)]; + +% if scheme == 2 % Ralston method +% +% b = [1/4 3/4]; +% +% k1 = f(told , uold); +% k2 = f(told + (2/3)*h, uold + (2/3)*h*(C\k1)); +% k = [transpose(k1); transpose(k2)]; +% kappa = [transpose(C\k1); transpose(C\k2)]; +% +% elseif scheme == 4 % 3/8th rule +% +% b = [1/8 3/8 3/8 1/8]; +% +% k1 = f(told , uold); +% k2 = f(told + (1/3)*h, uold + (1/3)*h*(C\k1)); +% k3 = f(told + (2/3)*h, uold - (1/3)*h*(C\k1) + h*(C\k2)); +% k4 = f(told + h , uold + h*(C\k1) - h*(C\k2) + h*(C\k3)); +% +% k = [transpose(k1); transpose(k2); transpose(k3); transpose(k4)]; +% +% elseif scheme == 8 % Dopri8 +% +% b = [(5.42937341165687622380535766363e-2) 0 0 0 0 (4.45031289275240888144113950566) (1.89151789931450038304281599044) -(5.8012039600105847814672114227) (3.1116436695781989440891606237e-1) -(1.52160949662516078556178806805e-1) (2.01365400804030348374776537501e-1) (4.47106157277725905176885569043e-2)]; +% +% k1 = f(told , uold); +% k2 = f(told + (0.526001519587677318785587544488e-01)*h , uold + (5.26001519587677318785587544488e-2)*h*(C\k1)); +% k3 = f(told + (0.789002279381515978178381316732e-01)*h , uold + (1.97250569845378994544595329183e-2)*h*(C\k1) + (5.91751709536136983633785987549e-2)*h*(C\k2)); +% k4 = f(told + (0.118350341907227396726757197510e+00)*h , uold + (2.95875854768068491816892993775e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (8.87627564304205475450678981324e-2)*h*(C\k3)); +% k5 = f(told + (0.281649658092772603273242802490e+00)*h , uold + (2.41365134159266685502369798665e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (-8.84549479328286085344864962717e-1)*h*(C\k3) + (9.24834003261792003115737966543e-1)*h*(C\k4)); +% k6 = f(told + (0.333333333333333333333333333333e+00)*h , uold + (3.7037037037037037037037037037e-2) *h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70828608729473871279604482173e-1)*h*(C\k4) + (1.25467687566822425016691814123e-1) *h*(C\k5)); +% k7 = f(told + (0.25e+00)*h , uold + (3.7109375e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70252211019544039314978060272e-1)*h*(C\k4) + (6.02165389804559606850219397283e-2) *h*(C\k5) + (-1.7578125e-2)*h*(C\k6)); +% k8 = f(told + (0.307692307692307692307692307692e+00)*h , uold + (3.70920001185047927108779319836e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70383925712239993810214054705e-1)*h*(C\k4) + (1.07262030446373284651809199168e-1) *h*(C\k5) + (-1.53194377486244017527936158236e-2)*h*(C\k6) + (8.27378916381402288758473766002e-3)*h*(C\k7)); +% k9 = f(told + (0.651282051282051282051282051282e+00)*h , uold + (6.24110958716075717114429577812e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-3.36089262944694129406857109825)*h*(C\k4) + (-8.68219346841726006818189891453e-1)*h*(C\k5) + (2.75920996994467083049415600797e1) *h*(C\k6) + (2.01540675504778934086186788979e1) *h*(C\k7) + (-4.34898841810699588477366255144e1)*h*(C\k8)); +% k10 = f(told + (0.6e+00)*h , uold + (4.77662536438264365890433908527e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-2.48811461997166764192642586468)*h*(C\k4) + (-5.90290826836842996371446475743e-1)*h*(C\k5) + (2.12300514481811942347288949897e1) *h*(C\k6) + (1.52792336328824235832596922938e1) *h*(C\k7) + (-3.32882109689848629194453265587e1)*h*(C\k8) + (-2.03312017085086261358222928593e-2) *h*(C\k9)); +% k11 = f(told + (0.857142857142857142857142857142e+00)*h, uold + (-9.3714243008598732571704021658e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (5.18637242884406370830023853209)*h*(C\k4) + (1.09143734899672957818500254654) *h*(C\k5) + (-8.14978701074692612513997267357) *h*(C\k6) + (-1.85200656599969598641566180701e1)*h*(C\k7) + (2.27394870993505042818970056734e1) *h*(C\k8) + (2.49360555267965238987089396762) *h*(C\k9) + (-3.0467644718982195003823669022) *h*(C\k10)); +% k12 = f(told + (0.1e+00)*h , uold + (2.27331014751653820792359768449) *h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-1.05344954667372501984066689879e1)*h*(C\k4) + (-2.00087205822486249909675718444) *h*(C\k5) + (-1.79589318631187989172765950534e1) *h*(C\k6) + (2.79488845294199600508499808837e1) *h*(C\k7) + (-2.85899827713502369474065508674) *h*(C\k8) + (-8.87285693353062954433549289258) *h*(C\k9) + (1.23605671757943030647266201528e1)*h*(C\k10) + (6.43392746015763530355970484046e-1)*h*(C\k11)); +% k = [transpose(k1); transpose(k2); transpose(k3); transpose(k4); transpose(k5); transpose(k6); transpose(k7); transpose(k8); transpose(k9); transpose(k10); transpose(k11); transpose(k12)]; +% else +% +% assert((scheme~=2)&&(scheme~=4)&&(scheme~=8) , 'You must pick a time scheme of 2nd, 4th or 8th order') +% +% end + + K = transpose(b*k); + K((n-nBCs+1):n) = 0; + unew = uold + C\(h*K); + tnew = told + h; + + j = j + 1; + u(:,j) = unew; + uold = unew; + t(j) = tnew; +end + +end diff --git a/source/@quantumGraph/timeEvolveARK.m b/source/@quantumGraph/timeEvolveARK.m new file mode 100644 index 00000000..53aa6950 --- /dev/null +++ b/source/@quantumGraph/timeEvolveARK.m @@ -0,0 +1,108 @@ +function [t,u] = timeEvolveARK(G,f,u0,h,tend) +% Given a graph G, function f and initial condition u0, this function +% evolves the solution to f(t,x) = u_xx in time using a 2nd, 4th or 8th +% order scheme until t=tend using the given step size h. The result will be +% a matrix the ith row of u is the solution evaluated at ith value of tvec. + +told = 0; % t_n +tnew = told; % t_n+1 +s = 3; % Number of b_i coeff's +j = 1; % Index for t and u +m = round(tend/h); % Number of t points +n = length(u0); % Number of x points + +nEdges = length(G.L); % Number of edges +nBCs = 2*nEdges; % Number of BC conditions +A = full(G.laplacianMatrix); % The second spatial derivative operator +B = full(G.weightMatrix); % The rectangular collocation or weight matrix +C = [B(1:(n-nBCs),:); A((n-nBCs+1):n,:)]; % Weight matrix diagonal with boundary data in bottom rows + + + +if size(u0,1)==1 % Corrects orientation of u0 + u0=u0'; +end + +u = zeros(n,m); +t = zeros(1,m); +u(:,1) = u0; +uold = u(:,1); % u_n + +x = 1.06858; +a = [x 0 0; 1/2-x x 0; 2*x 1-4*x x]; +b1 = 1/(6*(1-2*x)^2); +b2 = (3*(1-2*x)^2 - 1) / (3*(1-2*x)^2); +b3 = b1; +b = [b1 b2 b3]; +c = zeros(s,1); + +for i=1:s + c(i) = sum(a(i,:)); +end + + +while tnew < tend + told = tnew; + + k1guess = f(told , uold); + k1 = fsolve(@(x) x - f(told + c(1)*h, uold + a(1,1)*h*(C\x)), k1guess); + k2 = fsolve(@(x) x - f(told + c(2)*h, uold + a(2,1)*h*(C\k1) + a(2,2)*h*(C\x)), k1); + k3 = fsolve(@(x) x - f(told + c(3)*h, uold + a(3,1)*h*(C\k1) + a(3,2)*h*(C\k2) + a(3,3)*h*(C\x)), k2); + + k = [transpose(k1); transpose(k2); transpose(k3)]; + +% if scheme == 2 % Ralston method +% +% b = [1/4 3/4]; +% +% k1 = f(told , uold); +% k2 = f(told + (2/3)*h, uold + (2/3)*h*(C\k1)); +% k = [transpose(k1); transpose(k2)]; +% kappa = [transpose(C\k1); transpose(C\k2)]; +% +% elseif scheme == 4 % 3/8th rule +% +% b = [1/8 3/8 3/8 1/8]; +% +% k1 = f(told , uold); +% k2 = f(told + (1/3)*h, uold + (1/3)*h*(C\k1)); +% k3 = f(told + (2/3)*h, uold - (1/3)*h*(C\k1) + h*(C\k2)); +% k4 = f(told + h , uold + h*(C\k1) - h*(C\k2) + h*(C\k3)); +% +% k = [transpose(k1); transpose(k2); transpose(k3); transpose(k4)]; +% +% elseif scheme == 8 % Dopri8 +% +% b = [(5.42937341165687622380535766363e-2) 0 0 0 0 (4.45031289275240888144113950566) (1.89151789931450038304281599044) -(5.8012039600105847814672114227) (3.1116436695781989440891606237e-1) -(1.52160949662516078556178806805e-1) (2.01365400804030348374776537501e-1) (4.47106157277725905176885569043e-2)]; +% +% k1 = f(told , uold); +% k2 = f(told + (0.526001519587677318785587544488e-01)*h , uold + (5.26001519587677318785587544488e-2)*h*(C\k1)); +% k3 = f(told + (0.789002279381515978178381316732e-01)*h , uold + (1.97250569845378994544595329183e-2)*h*(C\k1) + (5.91751709536136983633785987549e-2)*h*(C\k2)); +% k4 = f(told + (0.118350341907227396726757197510e+00)*h , uold + (2.95875854768068491816892993775e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (8.87627564304205475450678981324e-2)*h*(C\k3)); +% k5 = f(told + (0.281649658092772603273242802490e+00)*h , uold + (2.41365134159266685502369798665e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (-8.84549479328286085344864962717e-1)*h*(C\k3) + (9.24834003261792003115737966543e-1)*h*(C\k4)); +% k6 = f(told + (0.333333333333333333333333333333e+00)*h , uold + (3.7037037037037037037037037037e-2) *h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70828608729473871279604482173e-1)*h*(C\k4) + (1.25467687566822425016691814123e-1) *h*(C\k5)); +% k7 = f(told + (0.25e+00)*h , uold + (3.7109375e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70252211019544039314978060272e-1)*h*(C\k4) + (6.02165389804559606850219397283e-2) *h*(C\k5) + (-1.7578125e-2)*h*(C\k6)); +% k8 = f(told + (0.307692307692307692307692307692e+00)*h , uold + (3.70920001185047927108779319836e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70383925712239993810214054705e-1)*h*(C\k4) + (1.07262030446373284651809199168e-1) *h*(C\k5) + (-1.53194377486244017527936158236e-2)*h*(C\k6) + (8.27378916381402288758473766002e-3)*h*(C\k7)); +% k9 = f(told + (0.651282051282051282051282051282e+00)*h , uold + (6.24110958716075717114429577812e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-3.36089262944694129406857109825)*h*(C\k4) + (-8.68219346841726006818189891453e-1)*h*(C\k5) + (2.75920996994467083049415600797e1) *h*(C\k6) + (2.01540675504778934086186788979e1) *h*(C\k7) + (-4.34898841810699588477366255144e1)*h*(C\k8)); +% k10 = f(told + (0.6e+00)*h , uold + (4.77662536438264365890433908527e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-2.48811461997166764192642586468)*h*(C\k4) + (-5.90290826836842996371446475743e-1)*h*(C\k5) + (2.12300514481811942347288949897e1) *h*(C\k6) + (1.52792336328824235832596922938e1) *h*(C\k7) + (-3.32882109689848629194453265587e1)*h*(C\k8) + (-2.03312017085086261358222928593e-2) *h*(C\k9)); +% k11 = f(told + (0.857142857142857142857142857142e+00)*h, uold + (-9.3714243008598732571704021658e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (5.18637242884406370830023853209)*h*(C\k4) + (1.09143734899672957818500254654) *h*(C\k5) + (-8.14978701074692612513997267357) *h*(C\k6) + (-1.85200656599969598641566180701e1)*h*(C\k7) + (2.27394870993505042818970056734e1) *h*(C\k8) + (2.49360555267965238987089396762) *h*(C\k9) + (-3.0467644718982195003823669022) *h*(C\k10)); +% k12 = f(told + (0.1e+00)*h , uold + (2.27331014751653820792359768449) *h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-1.05344954667372501984066689879e1)*h*(C\k4) + (-2.00087205822486249909675718444) *h*(C\k5) + (-1.79589318631187989172765950534e1) *h*(C\k6) + (2.79488845294199600508499808837e1) *h*(C\k7) + (-2.85899827713502369474065508674) *h*(C\k8) + (-8.87285693353062954433549289258) *h*(C\k9) + (1.23605671757943030647266201528e1)*h*(C\k10) + (6.43392746015763530355970484046e-1)*h*(C\k11)); +% k = [transpose(k1); transpose(k2); transpose(k3); transpose(k4); transpose(k5); transpose(k6); transpose(k7); transpose(k8); transpose(k9); transpose(k10); transpose(k11); transpose(k12)]; +% else +% +% assert((scheme~=2)&&(scheme~=4)&&(scheme~=8) , 'You must pick a time scheme of 2nd, 4th or 8th order') +% +% end + + K = transpose(b*k); + K((n-nBCs+1):n) = 0; + unew = uold + C\(h*K); + tnew = told + h; + + j = j + 1; + u(:,j) = unew; + uold = unew; + t(j) = tnew; +end + +end diff --git a/source/@quantumGraph/timeEvolveDAE.m b/source/@quantumGraph/timeEvolveDAE.m deleted file mode 100644 index ce1b2c81..00000000 --- a/source/@quantumGraph/timeEvolveDAE.m +++ /dev/null @@ -1,99 +0,0 @@ -function [u,t] = timeEvolveDAE(G,f,u0,h,tend,scheme) -% Given a graph G, function f and initial condition u0, this function -% evolves the solution to f(t,x) = u_xx in time using a 2nd, 4th or 8th -% order scheme until t=tend using the given step size h. The result will be -% a matrix the ith row of u is the solution evaluated at ith value of tvec. - -told = 0; % t_n -tnew = told; % t_n+1 -s = 4; % Number of b_i coeff's -j = 1; % Index for t and u -m = round(tend/h); % Number of t points -n = length(u0); % Number of x points - -nEdges = length(G.L); % Number of edges -nBCs = 2*nEdges; % Number of BC conditions -A = G.laplacianMatrix; % The second spatial derivative operator -B = G.weightMatrix; % The rectangular collocation or weight matrix -C = [B(1:(n-nBCs),:); A((n-nBCs+1):n,:)]; % Weight matrix diagonal with boundary data in bottom rows - -if size(u0,1)==1 % Corrects orientation of u0 - u0=u0'; -end - -u = zeros(n,m); -t = zeros(1,m); -u(:,1) = u0; -uold = u(:,1); % u_n -unew = uold; % u_n+1 - -a = [0 0 0 0; 1/2 0 0 0; 0 1/2 0 0; 0 0 1 0]; -c = zeros(s,1); - - -for i=1:s - c(i) = sum(a(i,:)); -end - - -while tnew < tend - told = tnew; - - if scheme == 2 % Ralston method - - b = [1/4 3/4]; - k = zeros(n,2); - - k1 = f(told , uold); - k2 = f(told + (2/3)*h, uold + (2/3)*h*(C\k1)); - k = [transpose(k1); transpose(k2)]; - - elseif scheme == 4 % 3/8th rule - - b = [1/8 3/8 3/8 1/8]; - k = zeros(n,4); - - k1 = f(told , uold); - k2 = f(told + (1/3)*h, uold + (1/2)*h*(C\k1)); - k3 = f(told + (2/3)*h, uold - (1/3)*h*(C\k1) + h*(C\k2)); - k4 = f(told + h , uold + h*(C\k1) - h*(C\k2) + h*(C\k3)); - k = [transpose(k1); transpose(k2); transpose(k3); transpose(k4)]; - - elseif scheme == 8 % Dopri8 - - b = [(5.42937341165687622380535766363e-2) 0 0 0 0 (4.45031289275240888144113950566) (1.89151789931450038304281599044) -(5.8012039600105847814672114227) (3.1116436695781989440891606237e-1) -(1.52160949662516078556178806805e-1) (2.01365400804030348374776537501e-1) (4.47106157277725905176885569043e-2)]; - k = zeros(n,12); - - k1 = f(told , uold); - k2 = f(told + (0.526001519587677318785587544488e-01)*h , uold + (5.26001519587677318785587544488e-2)*h*(C\k1)); - k3 = f(told + (0.789002279381515978178381316732e-01)*h , uold + (1.97250569845378994544595329183e-2)*h*(C\k1) + (5.91751709536136983633785987549e-2)*h*(C\k2)); - k4 = f(told + (0.118350341907227396726757197510e+00)*h , uold + (2.95875854768068491816892993775e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (8.87627564304205475450678981324e-2)*h*(C\k3)); - k5 = f(told + (0.281649658092772603273242802490e+00)*h , uold + (2.41365134159266685502369798665e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (-8.84549479328286085344864962717e-1)*h*(C\k3) + (9.24834003261792003115737966543e-1)*h*(C\k4)); - k6 = f(told + (0.333333333333333333333333333333e+00)*h , uold + (3.7037037037037037037037037037e-2) *h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70828608729473871279604482173e-1)*h*(C\k4) + (1.25467687566822425016691814123e-1) *h*(C\k5)); - k7 = f(told + (0.25e+00)*h , uold + (3.7109375e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70252211019544039314978060272e-1)*h*(C\k4) + (6.02165389804559606850219397283e-2) *h*(C\k5) + (-1.7578125e-2)*h*(C\k6)); - k8 = f(told + (0.307692307692307692307692307692e+00)*h , uold + (3.70920001185047927108779319836e-2)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (1.70383925712239993810214054705e-1)*h*(C\k4) + (1.07262030446373284651809199168e-1) *h*(C\k5) + (-1.53194377486244017527936158236e-2)*h*(C\k6) + (8.27378916381402288758473766002e-3)*h*(C\k7)); - k9 = f(told + (0.651282051282051282051282051282e+00)*h , uold + (6.24110958716075717114429577812e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-3.36089262944694129406857109825)*h*(C\k4) + (-8.68219346841726006818189891453e-1)*h*(C\k5) + (2.75920996994467083049415600797e1) *h*(C\k6) + (2.01540675504778934086186788979e1) *h*(C\k7) + (-4.34898841810699588477366255144e1)*h*(C\k8)); - k10 = f(told + (0.6e+00)*h , uold + (4.77662536438264365890433908527e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-2.48811461997166764192642586468)*h*(C\k4) + (-5.90290826836842996371446475743e-1)*h*(C\k5) + (2.12300514481811942347288949897e1) *h*(C\k6) + (1.52792336328824235832596922938e1) *h*(C\k7) + (-3.32882109689848629194453265587e1)*h*(C\k8) + (-2.03312017085086261358222928593e-2) *h*(C\k9)); - k11 = f(told + (0.857142857142857142857142857142e+00)*h, uold + (-9.3714243008598732571704021658e-1)*h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (5.18637242884406370830023853209)*h*(C\k4) + (1.09143734899672957818500254654) *h*(C\k5) + (-8.14978701074692612513997267357) *h*(C\k6) + (-1.85200656599969598641566180701e1)*h*(C\k7) + (2.27394870993505042818970056734e1) *h*(C\k8) + (2.49360555267965238987089396762) *h*(C\k9) + (-3.0467644718982195003823669022) *h*(C\k10)); - k12 = f(told + (0.1e+00)*h , uold + (2.27331014751653820792359768449) *h*(C\k1) + (0.0)*h*(C\k2) + (0.0)*h*(C\k3) + (-1.05344954667372501984066689879e1)*h*(C\k4) + (-2.00087205822486249909675718444) *h*(C\k5) + (-1.79589318631187989172765950534e1) *h*(C\k6) + (2.79488845294199600508499808837e1) *h*(C\k7) + (-2.85899827713502369474065508674) *h*(C\k8) + (-8.87285693353062954433549289258) *h*(C\k9) + (1.23605671757943030647266201528e1)*h*(C\k10) + (6.43392746015763530355970484046e-1)*h*(C\k11)); - k = [transpose(k1); transpose(k2); transpose(k3); transpose(k4); transpose(k5); transpose(k6); transpose(k7); transpose(k8); transpose(k9); transpose(k10); transpose(k11); transpose(k12)]; - - else - - error('You must pick a time scheme of 2nd, 4th or 8th order') - - end - - K = transpose(b*k); - unew = C \ (B*uold + h*K); - tnew = told + h; - - j = j + 1; - u(:,j) = unew; - uold = unew; - t(j) = tnew; -end - - - -end \ No newline at end of file diff --git a/source/chebyshev/chebyshevDeriv.m b/source/chebyshev/chebyshevDeriv.m index c4f839f8..f96f1aa9 100644 --- a/source/chebyshev/chebyshevDeriv.m +++ b/source/chebyshev/chebyshevDeriv.m @@ -15,5 +15,5 @@ dX = X-X'; D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries D = D - diag(sum(D')); % diagonal entries -D = 2/l*D; % rescaling for length of interval +D = 2/l*D; % rescale for length of interval end \ No newline at end of file diff --git a/source/chebyshev/clencurt.m b/source/chebyshev/clencurt.m index 9f413bb3..aad17a05 100644 --- a/source/chebyshev/clencurt.m +++ b/source/chebyshev/clencurt.m @@ -5,22 +5,14 @@ % with a Chebyshev discretization. Also produces the weights and % coefficients used in the Clenshaw-Curtis quadrature. -n = length(f); % Total number of discretization points (including end points?) +n = length(f); % Total number of discretization points including end points N = n-1; - -% % WEIGHTS: fast version I DONT THINK I NEED THE WEIGHTS!!!!!! -% c = zeros(n,2); -% c(1:2:n,1) = (2./[1 1-(2:2:N).^2 ])'; c(2,2) = 1; -% f1 = real(ifft([c(1:n,:);c(N:-1:2,:)])); -% w = L*([f1(1,1); 2*f1(2:N,1); f1(n,1)])/2; - - % COEFFICENTS: a_k = int_0^pi f(cos(x))*cos(kx) dx theta = pi*(0:N)/N; dtheta = theta(2) - theta(1); -a0 = 2/pi * (dot(f,cos(0*theta)) - .5 *(f(1)*cos(0*theta(1)) + f(n)*cos(0*theta(n))) ) * dtheta; +a0 = 2/pi * (dot(f,cos(0*theta)) - .5 *(f(1)*cos(0*theta(1)) + f(n)*cos(0*theta(n)))) * dtheta; z = a0; eps = 10^(-16); maxk = 10^4; @@ -28,7 +20,7 @@ for k=2:maxk if mod(k,2)==0 % Picks out even indicies - ak = 2/pi * (dot(f,cos(k*theta)) - .5 *(f(1)*cos(k*theta(1)) + f(n)*cos(k*theta(n))) ) * dtheta; % Trapezoid rule + ak = 2/pi * (dot(f,cos(k*theta)) - .5 *(f(1)*cos(k*theta(1)) + f(n)*cos(k*theta(n)))) * dtheta; % Trapezoid rule bk = 2*ak/(1-k^2); if abs(bk)>eps && k N2 + [PhiBif,muBif] = HHBLocator(funcs, [PhiCol1';PhiCol2'], [mu1;mu2]); + else + assert(~isnan(PhiBif),'Failed to find any HHB.') + end +end diff --git a/source/continuation/acm/HHBLocator.m b/source/continuation/acm/HHBLocator.m new file mode 100644 index 00000000..2b541548 --- /dev/null +++ b/source/continuation/acm/HHBLocator.m @@ -0,0 +1,75 @@ +function [PhiBif,muBif] = HHBLocator(fcns, PhiColVec, muVec) +% After determing that there is a bifurcation point between muVec[1] and +% muVec[2], this program is then used to locate the bifurcation point muBif +% and determine its affiliated eigenfunction, PhiBif. +% +% Parameters +% ---------- +% fcns : all necessary functions associated with the cubic NLS +% muVec : initial eigenvalues on either side of HHB +% PhiColVec : stationary state solutions associated with muVec +% +% Outputs +% ------- +% PhiBif : the NLS solution at the HHB +% muBif : the HHB point + + mu1 = muVec(1); + mu2 = muVec(2); + PhiCol1 = PhiColVec(:,1); + PhiCol2 = PhiColVec(:,2); + muGuess = (mu1 + mu2)/2; % Initial guess for HHB point + muBif = muGuess; + error = 10^(-6); + + myF = @(z,mu) fcns.f(z,mu); + myMatrix = @(u) fcns.fLinMatrix(u,mu); + if abs(muGuess-mu1) < abs(mu2-muGuess) % Checks to see if muGuess is closer to mu1 than mu2 + myFGuess = @(z) myF(z,muGuess); + myMatrixGuess = @(u) myMatrix(u,muGuess); + [PhiGuess,~,~] = solveNewton(PhiCol1,myFGuess,myMatrixGuess,); + else + myFGuess = @(z) myF(z,muGuess); + myMatrixGuess = @(u) myMatrix(u,muGuess); + [PhiGuess,~,~] = solveNewton(PhiCol2,myFGuess,myMatrixGuess,error); + end + PhiBif = PhiGuess; % Initial guess for solution at muGuess + + residual = abs(mu1-mu2); + error = 1e-8; + maxTries = 200; + tries = 0; + + while residual > error && tries < maxTries + tries = tries+1; + + if HHBDetector(Phi,[PhiCol1,PhiGuess],[mu1,muGuess]) == 1 % When the bifurcation is between mu1 and muGuess + mu2 = muGuess; + PhiCol2 = PhiGuess; + else % When the bifurcation is between muGuess and mu2 + mu1 = muGuess; + PhiCol1 = PhiGuess; + end + + muGuess = (mu1 + mu2)/2; + + myF = @(z,mu) fcns.f(z,mu); + myMatrix = @(u) fcns.fLinMatrix(u,mu); + if abs(muGuess-mu1) < abs(m2-muGuess) + myFGuess = @(z) myF(z,muGuess); + myMatrixGuess = @(u) myMatrix(u,muGuess); + [PhiGuess,~,~] = solveNewton(PhiCol1,myFGuess,myMatrixGuess,error); + else + myFGuess = @(z) myF(z,muGuess); + myMatrixGuess = @(u) myMatrix(u,muGuess); + [PhiGuess,~,~] = solveNewton(PhiCol2,myFGuess,myMatrixGuess,error); + end + + residual = abs(mu1-mu2); + muBif = muGuess; + PhiBif = PhiGuess; + end + + assert(triesNormal*) (* Beginning of Notebook Content *) @@ -56,8 +56,7 @@ Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{ - SubscriptBox["y", "2"], "''"}], "[", "x", "]"}], "==", "x"}]}]], "Input",\ - + SubscriptBox["y", "2"], "''"}], "[", "x", "]"}], "==", "x"}]}]], "Input", CellChangeTimes->{{3.834678721172653*^9, 3.8346787390175247`*^9}}, CellLabel->"In[2]:=",ExpressionUUID->"af038ae0-03e8-4ca0-b062-9786503f1d09"], @@ -103,8 +102,7 @@ Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{ - SubscriptBox["y", "4"], "''"}], "[", "x", "]"}], "==", "1"}]}]], "Input",\ - + SubscriptBox["y", "4"], "''"}], "[", "x", "]"}], "==", "1"}]}]], "Input", CellChangeTimes->{{3.8346787573015833`*^9, 3.8346787806882544`*^9}}, CellLabel->"In[4]:=",ExpressionUUID->"0d88c929-92a4-462b-bb39-1be0958481dd"], @@ -137,7 +135,7 @@ Cell[BoxData[ RowBox[{ SubscriptBox["y", "2"], "'"}], "[", "0", "]"}], "+", RowBox[{ - SubscriptBox["y", "1"], "[", "0", "]"}]}], "==", "1"}], ",", + SubscriptBox["y", "1"], "[", "0", "]"}]}], "\[Equal]", "1"}], ",", RowBox[{ RowBox[{ SubscriptBox["y", "1"], "[", "0", "]"}], "==", @@ -147,8 +145,10 @@ Cell[BoxData[ RowBox[{ SubscriptBox["y", "2"], "[", "0", "]"}]}]}], "}"}]}]], "Input", CellChangeTimes->{{3.834678815307212*^9, 3.834678960354398*^9}, - 3.834679091815206*^9, {3.834680101997093*^9, 3.834680108768464*^9}}, - CellLabel->"In[10]:=",ExpressionUUID->"8e1435ad-73a8-410c-a837-1ec335ea5781"], + 3.834679091815206*^9, {3.834680101997093*^9, 3.834680108768464*^9}, { + 3.8506535281132536`*^9, 3.850653528401059*^9}, + 3.850654240775215*^9},ExpressionUUID->"8e1435ad-73a8-410c-a837-\ +1ec335ea5781"], Cell[BoxData[ RowBox[{"{", @@ -205,7 +205,7 @@ Cell[BoxData[ RowBox[{ SubscriptBox["y", "4"], "'"}], "[", "0", "]"}], "+", RowBox[{ - SubscriptBox["y", "3"], "[", "0", "]"}]}], "==", "2"}], ",", + SubscriptBox["y", "3"], "[", "0", "]"}]}], "\[Equal]", "2"}], ",", RowBox[{ RowBox[{ SubscriptBox["y", "2"], "[", "4", "]"}], "==", @@ -217,9 +217,10 @@ Cell[BoxData[ RowBox[{ SubscriptBox["y", "4"], "[", "0", "]"}]}]}], "}"}]}]], "Input", CellChangeTimes->{{3.834678964818534*^9, 3.834679054548163*^9}, { - 3.83467909814336*^9, 3.8346791212260847`*^9}, {3.834680117471212*^9, - 3.834680131890923*^9}}, - CellLabel->"In[11]:=",ExpressionUUID->"6048eb9b-b2ac-43b3-8d15-e297dddfe78d"], + 3.83467909814336*^9, 3.8346791212260847`*^9}, {3.834680117471212*^9, + 3.834680131890923*^9}, 3.8506535259593425`*^9, + 3.8506542366069365`*^9},ExpressionUUID->"6048eb9b-b2ac-43b3-8d15-\ +e297dddfe78d"], Cell[BoxData[ RowBox[{"{", @@ -265,9 +266,11 @@ Cell[BoxData[ RowBox[{"Vertex3", "=", " ", RowBox[{ RowBox[{ - SubscriptBox["y", "4"], "[", "1", "]"}], "==", "3"}]}]], "Input", - CellChangeTimes->{{3.834679129823811*^9, 3.8346791514351797`*^9}}, - CellLabel->"In[8]:=",ExpressionUUID->"f855c0ba-8ebe-4922-aea4-168ec05af293"], + SubscriptBox["y", "4"], "[", "1", "]"}], "\[Equal]", "3"}]}]], "Input", + CellChangeTimes->{{3.834679129823811*^9, 3.8346791514351797`*^9}, { + 3.8506535315219784`*^9, 3.8506535323017616`*^9}, + 3.850654230841838*^9},ExpressionUUID->"f855c0ba-8ebe-4922-aea4-\ +168ec05af293"], Cell[BoxData[ RowBox[{ @@ -352,9 +355,9 @@ Cell[BoxData[ CellLabel->"Out[12]=",ExpressionUUID->"955ef5bf-cb5b-4c11-863b-bb40aa02e28c"] }, Open ]] }, -WindowSize->{808, 747}, -WindowMargins->{{Automatic, -1108}, {Automatic, -267}}, -FrontEndVersion->"12.3 for Mac OS X x86 (64-bit) (May 11, 2021)", +WindowSize->{1026., 622.5}, +WindowMargins->{{-4.875, Automatic}, {Automatic, -4.875}}, +FrontEndVersion->"12.1 for Microsoft Windows (64-bit) (June 19, 2020)", StyleDefinitions->"Default.nb", ExpressionUUID->"526e867f-66c1-422b-ac1c-bcab39534acc" ] @@ -371,36 +374,36 @@ CellTagsIndex->{} Notebook[{ Cell[558, 20, 241, 6, 35, "Text",ExpressionUUID->"66afe82f-2a44-4e28-8366-15904e016d78"], Cell[CellGroupData[{ -Cell[824, 30, 327, 8, 30, "Input",ExpressionUUID->"3434e81c-0d26-41f2-a47c-2add7e5727f7"], -Cell[1154, 40, 328, 8, 34, "Output",ExpressionUUID->"bc475ffd-bb1f-4c6e-935b-a074dbf16bee"] +Cell[824, 30, 327, 8, 28, "Input",ExpressionUUID->"3434e81c-0d26-41f2-a47c-2add7e5727f7"], +Cell[1154, 40, 328, 8, 32, "Output",ExpressionUUID->"bc475ffd-bb1f-4c6e-935b-a074dbf16bee"] }, Open ]], Cell[CellGroupData[{ -Cell[1519, 53, 300, 8, 30, "Input",ExpressionUUID->"af038ae0-03e8-4ca0-b062-9786503f1d09"], -Cell[1822, 63, 300, 7, 34, "Output",ExpressionUUID->"d7c873ad-6983-43f6-81aa-dcafb2bf389c"] +Cell[1519, 53, 298, 7, 28, "Input",ExpressionUUID->"af038ae0-03e8-4ca0-b062-9786503f1d09"], +Cell[1820, 62, 300, 7, 32, "Output",ExpressionUUID->"d7c873ad-6983-43f6-81aa-dcafb2bf389c"] }, Open ]], Cell[CellGroupData[{ -Cell[2159, 75, 347, 9, 30, "Input",ExpressionUUID->"11d1673b-7c38-442e-8215-2e7acba2ca70"], -Cell[2509, 86, 352, 9, 34, "Output",ExpressionUUID->"663d05aa-a2c7-459f-b28a-f8aa289f7461"] +Cell[2157, 74, 347, 9, 28, "Input",ExpressionUUID->"11d1673b-7c38-442e-8215-2e7acba2ca70"], +Cell[2507, 85, 352, 9, 32, "Output",ExpressionUUID->"663d05aa-a2c7-459f-b28a-f8aa289f7461"] }, Open ]], Cell[CellGroupData[{ -Cell[2898, 100, 302, 8, 30, "Input",ExpressionUUID->"0d88c929-92a4-462b-bb39-1be0958481dd"], -Cell[3203, 110, 298, 7, 34, "Output",ExpressionUUID->"913c40fe-a765-4bd6-bd5f-7d7d14087e6a"] +Cell[2896, 99, 300, 7, 28, "Input",ExpressionUUID->"0d88c929-92a4-462b-bb39-1be0958481dd"], +Cell[3199, 108, 298, 7, 32, "Output",ExpressionUUID->"913c40fe-a765-4bd6-bd5f-7d7d14087e6a"] }, Open ]], Cell[CellGroupData[{ -Cell[3538, 122, 960, 28, 30, "Input",ExpressionUUID->"8e1435ad-73a8-410c-a837-1ec335ea5781"], -Cell[4501, 152, 1000, 29, 34, "Output",ExpressionUUID->"bbeec3fa-6651-4426-b9bc-16807fc03791"] +Cell[3534, 120, 1022, 30, 28, "Input",ExpressionUUID->"8e1435ad-73a8-410c-a837-1ec335ea5781"], +Cell[4559, 152, 1000, 29, 32, "Output",ExpressionUUID->"bbeec3fa-6651-4426-b9bc-16807fc03791"] }, Open ]], Cell[CellGroupData[{ -Cell[5538, 186, 1173, 35, 30, "Input",ExpressionUUID->"6048eb9b-b2ac-43b3-8d15-e297dddfe78d"], -Cell[6714, 223, 1204, 35, 34, "Output",ExpressionUUID->"054bbe80-9494-495e-9507-55959b032760"] +Cell[5596, 186, 1211, 36, 28, "Input",ExpressionUUID->"6048eb9b-b2ac-43b3-8d15-e297dddfe78d"], +Cell[6810, 224, 1204, 35, 32, "Output",ExpressionUUID->"054bbe80-9494-495e-9507-55959b032760"] }, Open ]], Cell[CellGroupData[{ -Cell[7955, 263, 283, 6, 30, "Input",ExpressionUUID->"f855c0ba-8ebe-4922-aea4-168ec05af293"], -Cell[8241, 271, 227, 5, 34, "Output",ExpressionUUID->"e1a3f17f-d98f-468b-97f2-375f8e5f14e5"] +Cell[8051, 264, 348, 8, 28, "Input",ExpressionUUID->"f855c0ba-8ebe-4922-aea4-168ec05af293"], +Cell[8402, 274, 227, 5, 32, "Output",ExpressionUUID->"e1a3f17f-d98f-468b-97f2-375f8e5f14e5"] }, Open ]], Cell[CellGroupData[{ -Cell[8505, 281, 550, 14, 52, "Input",ExpressionUUID->"521264de-eb40-45e0-8e82-a8919d380f73"], -Cell[9058, 297, 1730, 54, 90, "Output",ExpressionUUID->"955ef5bf-cb5b-4c11-863b-bb40aa02e28c"] +Cell[8666, 284, 550, 14, 28, "Input",ExpressionUUID->"521264de-eb40-45e0-8e82-a8919d380f73"], +Cell[9219, 300, 1730, 54, 85, "Output",ExpressionUUID->"955ef5bf-cb5b-4c11-863b-bb40aa02e28c"] }, Open ]] } ] diff --git a/source/examples/poissonExampleDumbbellNonhomogeneous3.mlx b/source/examples/poissonExampleDumbbellNonhomogeneous3.mlx new file mode 100644 index 00000000..156e6f9a Binary files /dev/null and b/source/examples/poissonExampleDumbbellNonhomogeneous3.mlx differ diff --git a/source/examples/poissonExampleLollipop.nb b/source/examples/poissonExampleLollipop.nb new file mode 100644 index 00000000..273c3015 --- /dev/null +++ b/source/examples/poissonExampleLollipop.nb @@ -0,0 +1,631 @@ +(* Content-type: application/vnd.wolfram.mathematica *) + +(*** Wolfram Notebook File ***) +(* http://www.wolfram.com/nb *) + +(* CreatedBy='Mathematica 12.1' *) + +(*CacheID: 234*) +(* Internal cache information: +NotebookFileLineBreakTest +NotebookFileLineBreakTest +NotebookDataPosition[ 158, 7] +NotebookDataLength[ 23235, 621] +NotebookOptionsPosition[ 20386, 561] +NotebookOutlinePosition[ 20788, 577] +CellTagsIndexPosition[ 20745, 574] +WindowFrame->Normal*) + +(* Beginning of Notebook Content *) +Notebook[{ + +Cell[CellGroupData[{ +Cell[TextData[{ + "Solve a Poisson problem on a quantum graph subject to ", + StyleBox["homogeneous", + FontWeight->"Bold"], + " vertex conditions." +}], "Subsection", + CellChangeTimes->{{3.834758025487357*^9, 3.8347580710539103`*^9}, + 3.850657970006077*^9},ExpressionUUID->"459730e3-b4c0-4dc3-86d6-\ +d1c4affa7b98"], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Side1", "=", + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "''"}], "[", "x", "]"}], "\[Equal]", + RowBox[{ + RowBox[{"-", + SuperscriptBox["\[Pi]", "2"]}], + RowBox[{"Cos", "[", + RowBox[{"\[Pi]", "*", "x"}], "]"}]}]}]}]], "Input", + CellChangeTimes->{{3.834678685914446*^9, 3.834678717948956*^9}, { + 3.8506567700859294`*^9, 3.8506567710683413`*^9}, {3.850657077840087*^9, + 3.850657086836153*^9}, 3.8506571199286633`*^9, {3.850657227765092*^9, + 3.850657238043795*^9}, {3.850657362094737*^9, 3.8506573676940765`*^9}}, + CellLabel-> + "In[119]:=",ExpressionUUID->"c2ce4243-f67f-433e-af2d-d14d5a6e465c"], + +Cell[BoxData[ + RowBox[{ + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "1"], "\[Prime]\[Prime]", + MultilineFunction->None], "[", "x", "]"}], "\[Equal]", + RowBox[{ + RowBox[{"-", + SuperscriptBox["\[Pi]", "2"]}], " ", + RowBox[{"Cos", "[", + RowBox[{"\[Pi]", " ", "x"}], "]"}]}]}]], "Output", + CellChangeTimes->{ + 3.834678720567547*^9, 3.850654246431823*^9, 3.8506546577010136`*^9, + 3.850655599120631*^9, 3.8506557591957493`*^9, 3.850656201705437*^9, + 3.850656308776461*^9, 3.8506564236377*^9, 3.850656579005994*^9, + 3.8506567052850513`*^9, {3.8506567741580057`*^9, 3.850656775645717*^9}, + 3.850656892847095*^9, 3.850656949453047*^9, 3.850657008693903*^9, { + 3.85065709523395*^9, 3.850657124680691*^9}, 3.8506573427306414`*^9, { + 3.85065738002584*^9, 3.8506574148846154`*^9}, 3.8506581159328227`*^9}, + CellLabel-> + "Out[119]=",ExpressionUUID->"a0571591-7e57-47a7-8111-76c01f353f9a"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Side2", "=", + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "''"}], "[", "x", "]"}], "\[Equal]", + RowBox[{ + RowBox[{"-", "4"}], + RowBox[{"Cos", "[", + RowBox[{"2", "x"}], "]"}]}]}]}]], "Input", + CellChangeTimes->{{3.834678721172653*^9, 3.8346787390175247`*^9}, { + 3.850656062538176*^9, 3.8506560643083076`*^9}, {3.8506566999446516`*^9, + 3.8506567007639923`*^9}, 3.850657001866477*^9, 3.8506570912282743`*^9, { + 3.8506572298961353`*^9, 3.850657241322625*^9}, {3.8506573695540805`*^9, + 3.8506573757934012`*^9}, 3.850657411925483*^9}, + CellLabel-> + "In[120]:=",ExpressionUUID->"b1e599e2-c453-4afb-b14b-eb0509d45b0e"], + +Cell[BoxData[ + RowBox[{ + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "2"], "\[Prime]\[Prime]", + MultilineFunction->None], "[", "x", "]"}], "\[Equal]", + RowBox[{ + RowBox[{"-", "4"}], " ", + RowBox[{"Cos", "[", + RowBox[{"2", " ", "x"}], "]"}]}]}]], "Output", + CellChangeTimes->{ + 3.850657124786903*^9, 3.8506573427936964`*^9, {3.8506573800984526`*^9, + 3.8506574149535217`*^9}, 3.850658116024008*^9}, + CellLabel-> + "Out[120]=",ExpressionUUID->"b1330e79-85cb-4ddf-9406-88895e4a7d40"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Vertex1", "=", + RowBox[{"{", + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "'"}], "[", "0", "]"}], "\[Equal]", "0"}], + "}"}]}]], "Input", + CellChangeTimes->{{3.834678964818534*^9, 3.834679054548163*^9}, { + 3.83467909814336*^9, 3.8346791212260847`*^9}, {3.834680117471212*^9, + 3.834680131890923*^9}, 3.8506535259593425`*^9, 3.8506542366069365`*^9, + 3.850654566481948*^9, {3.850654609239808*^9, 3.8506546319391756`*^9}, { + 3.850655754862317*^9, 3.850655755260985*^9}, {3.8506561282878056`*^9, + 3.850656128506816*^9}, {3.850656179326333*^9, 3.850656194903241*^9}, { + 3.8506562948960257`*^9, 3.850656298173747*^9}, {3.8506568182908597`*^9, + 3.8506568253680487`*^9}, {3.85065687610748*^9, 3.8506568891480427`*^9}, { + 3.850657251237673*^9, 3.8506572879239845`*^9}}, + CellLabel-> + "In[121]:=",ExpressionUUID->"2ca5a43f-0e7c-4d8f-a01a-3d97f0f7c774"], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "1"], "\[Prime]", + MultilineFunction->None], "[", "0", "]"}], "\[Equal]", "0"}], + "}"}]], "Output", + CellChangeTimes->{ + 3.85065642381654*^9, 3.85065657920348*^9, 3.8506567054719696`*^9, + 3.8506567758191576`*^9, 3.8506568930209246`*^9, 3.8506569497505836`*^9, + 3.850657008890003*^9, {3.85065709545033*^9, 3.8506571249798975`*^9}, + 3.850657342905113*^9, {3.8506573802064686`*^9, 3.8506574150084963`*^9}, + 3.850658116093482*^9}, + CellLabel-> + "Out[121]=",ExpressionUUID->"a697f319-b098-49f5-afe1-f7eb7fc72992"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Vertex2", "=", + RowBox[{"{", + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "'"}], "[", "0", "]"}], "-", + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "'"}], "[", "\[Pi]", "]"}], "+", + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "'"}], "[", "2", "]"}]}], "\[Equal]", "0"}], + ",", + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "[", "0", "]"}], "==", + RowBox[{ + SubscriptBox["y", "2"], "[", "\[Pi]", "]"}], "==", + RowBox[{ + SubscriptBox["y", "1"], "[", "2", "]"}]}]}], "}"}]}]], "Input", + CellChangeTimes->{{3.834678815307212*^9, 3.834678960354398*^9}, + 3.834679091815206*^9, {3.834680101997093*^9, 3.834680108768464*^9}, { + 3.8506535281132536`*^9, 3.850653528401059*^9}, 3.850654240775215*^9, { + 3.8506546241912746`*^9, 3.850654625032473*^9}, {3.8506555832915897`*^9, + 3.850655587715314*^9}, {3.850656105505641*^9, 3.8506561059986706`*^9}, + 3.850656150647132*^9, {3.8506562839373245`*^9, 3.850656287632543*^9}, { + 3.850656791161252*^9, 3.850656812949095*^9}, {3.850657249495758*^9, + 3.8506572496053743`*^9}, {3.850657299623186*^9, 3.850657336922967*^9}}, + CellLabel-> + "In[122]:=",ExpressionUUID->"fbbf2e80-0c95-4928-bfff-f47af7a9fe52"], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "1"], "\[Prime]", + MultilineFunction->None], "[", "2", "]"}], "+", + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "2"], "\[Prime]", + MultilineFunction->None], "[", "0", "]"}], "-", + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "2"], "\[Prime]", + MultilineFunction->None], "[", "\[Pi]", "]"}]}], "\[Equal]", "0"}], + ",", + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "[", "0", "]"}], "\[Equal]", + RowBox[{ + SubscriptBox["y", "2"], "[", "\[Pi]", "]"}], "\[Equal]", + RowBox[{ + SubscriptBox["y", "1"], "[", "2", "]"}]}]}], "}"}]], "Output", + CellChangeTimes->{ + 3.8346789636705008`*^9, 3.834680109930167*^9, 3.85065424687862*^9, + 3.8506546578932705`*^9, 3.850655599313405*^9, 3.850655759382452*^9, + 3.8506562018607187`*^9, 3.8506563089355545`*^9, 3.850656423759699*^9, + 3.850656579146865*^9, 3.850656705414361*^9, 3.8506567757631865`*^9, + 3.850656892961679*^9, 3.8506569496581817`*^9, 3.8506570088279743`*^9, { + 3.8506570953848543`*^9, 3.85065712492*^9}, {3.8506573411184406`*^9, + 3.8506573429656816`*^9}, {3.850657380269035*^9, 3.8506574150704403`*^9}, + 3.850658116155764*^9}, + CellLabel-> + "Out[122]=",ExpressionUUID->"782125e5-9d48-4f0c-85a8-170051501035"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"DSolve", "[", + RowBox[{ + RowBox[{"{", + RowBox[{"Side1", ",", "Side2", ",", "Vertex1", ",", "Vertex2"}], "}"}], + ",", + RowBox[{"{", + RowBox[{ + SubscriptBox["y", "1"], ",", + SubscriptBox["y", "2"]}], "}"}], ",", "x"}], "]"}]], "Input", + CellChangeTimes->{{3.834679159495409*^9, 3.834679199735717*^9}, { + 3.85065464029145*^9, 3.850654651576705*^9}}, + CellLabel-> + "In[123]:=",ExpressionUUID->"f6d4efcc-c2fd-430f-a82d-c896e07895ac"], + +Cell[BoxData[ + TemplateBox[{ + "DSolve", "bvsing", + "\"Unable to resolve some of the arbitrary constants in the general \ +solution using the given boundary conditions. It is possible that some of the \ +conditions have been specified at a singular point for the equation.\"", 2, + 123, 15, 28936614898608263030, "Local"}, + "MessageTemplate"]], "Message", "MSG", + CellChangeTimes->{ + 3.8506563090696154`*^9, 3.850656423940381*^9, 3.85065657933856*^9, + 3.850656705611519*^9, 3.850656775957872*^9, 3.8506568930700264`*^9, + 3.850656949842327*^9, 3.8506570089502172`*^9, {3.850657095521228*^9, + 3.8506571250387325`*^9}, 3.850657343033668*^9, {3.850657380333579*^9, + 3.850657415137489*^9}, 3.8506581162314157`*^9}, + CellLabel-> + "During evaluation of \ +In[123]:=",ExpressionUUID->"fba1ea3b-1021-40d6-8754-0dd930352e64"], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{"{", + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "\[Rule]", + RowBox[{"Function", "[", + RowBox[{ + RowBox[{"{", "x", "}"}], ",", + RowBox[{ + TemplateBox[{"1"}, + "C"], "+", + RowBox[{"Cos", "[", + RowBox[{"\[Pi]", " ", "x"}], "]"}]}]}], "]"}]}], ",", + RowBox[{ + SubscriptBox["y", "2"], "\[Rule]", + RowBox[{"Function", "[", + RowBox[{ + RowBox[{"{", "x", "}"}], ",", + RowBox[{ + TemplateBox[{"1"}, + "C"], "+", + RowBox[{"Cos", "[", + RowBox[{"2", " ", "x"}], "]"}]}]}], "]"}]}]}], "}"}], + "}"}]], "Output", + CellChangeTimes->{ + 3.8346792010932693`*^9, 3.8346801435371943`*^9, 3.8506542471161647`*^9, + 3.8506546580415635`*^9, 3.8506555994399376`*^9, 3.850655759509415*^9, + 3.8506562019967995`*^9, 3.850656309179753*^9, 3.850656423951334*^9, + 3.85065657934756*^9, 3.8506567056216993`*^9, 3.8506567759728107`*^9, + 3.8506568930810633`*^9, 3.8506569498532786`*^9, 3.8506570089621677`*^9, { + 3.850657095533163*^9, 3.850657125050787*^9}, 3.850657343041707*^9, { + 3.8506573803445215`*^9, 3.8506574151474514`*^9}, 3.850658116244913*^9}, + CellLabel-> + "Out[123]=",ExpressionUUID->"7dc5a94a-4de5-4e7c-be83-fcc6e20b3b3d"] +}, Open ]] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[TextData[{ + "Solve a Poisson problem on a quantum graph subject to ", + StyleBox["nonhomogeneous", + FontWeight->"Bold"], + " vertex conditions." +}], "Subsection", + CellChangeTimes->{{3.834758025487357*^9, 3.8347580710539103`*^9}, + 3.850657970006077*^9, {3.8506580030622573`*^9, + 3.850658007012293*^9}},ExpressionUUID->"6ffa0d99-169a-43b4-89fe-\ +f3d442d8d5fc"], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Side1", "=", + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "''"}], "[", "x", "]"}], "\[Equal]", + RowBox[{ + RowBox[{"-", + SuperscriptBox["\[Pi]", "2"]}], + RowBox[{"Cos", "[", + RowBox[{"\[Pi]", "*", "x"}], "]"}]}]}]}]], "Input", + CellChangeTimes->{{3.834678685914446*^9, 3.834678717948956*^9}, { + 3.8506567700859294`*^9, 3.8506567710683413`*^9}, {3.850657077840087*^9, + 3.850657086836153*^9}, 3.8506571199286633`*^9, {3.850657227765092*^9, + 3.850657238043795*^9}, {3.850657362094737*^9, 3.8506573676940765`*^9}}, + CellLabel-> + "In[134]:=",ExpressionUUID->"a25ca004-781d-4237-bb37-d5989eb275ed"], + +Cell[BoxData[ + RowBox[{ + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "1"], "\[Prime]\[Prime]", + MultilineFunction->None], "[", "x", "]"}], "\[Equal]", + RowBox[{ + RowBox[{"-", + SuperscriptBox["\[Pi]", "2"]}], " ", + RowBox[{"Cos", "[", + RowBox[{"\[Pi]", " ", "x"}], "]"}]}]}]], "Output", + CellChangeTimes->{ + 3.834678720567547*^9, 3.850654246431823*^9, 3.8506546577010136`*^9, + 3.850655599120631*^9, 3.8506557591957493`*^9, 3.850656201705437*^9, + 3.850656308776461*^9, 3.8506564236377*^9, 3.850656579005994*^9, + 3.8506567052850513`*^9, {3.8506567741580057`*^9, 3.850656775645717*^9}, + 3.850656892847095*^9, 3.850656949453047*^9, 3.850657008693903*^9, { + 3.85065709523395*^9, 3.850657124680691*^9}, 3.8506573427306414`*^9, { + 3.85065738002584*^9, 3.8506574148846154`*^9}, 3.850658116314929*^9, + 3.850658150934472*^9, 3.8506581913655224`*^9}, + CellLabel-> + "Out[134]=",ExpressionUUID->"41e24e3d-6bb6-4743-b3b4-e1702e2bbd7a"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Side2", "=", + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "''"}], "[", "x", "]"}], "\[Equal]", + RowBox[{ + RowBox[{"-", "4"}], + RowBox[{"Cos", "[", + RowBox[{"2", "x"}], "]"}]}]}]}]], "Input", + CellChangeTimes->{{3.834678721172653*^9, 3.8346787390175247`*^9}, { + 3.850656062538176*^9, 3.8506560643083076`*^9}, {3.8506566999446516`*^9, + 3.8506567007639923`*^9}, 3.850657001866477*^9, 3.8506570912282743`*^9, { + 3.8506572298961353`*^9, 3.850657241322625*^9}, {3.8506573695540805`*^9, + 3.8506573757934012`*^9}, 3.850657411925483*^9}, + CellLabel-> + "In[135]:=",ExpressionUUID->"1a990e1b-45f5-4a75-b891-c354571430b0"], + +Cell[BoxData[ + RowBox[{ + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "2"], "\[Prime]\[Prime]", + MultilineFunction->None], "[", "x", "]"}], "\[Equal]", + RowBox[{ + RowBox[{"-", "4"}], " ", + RowBox[{"Cos", "[", + RowBox[{"2", " ", "x"}], "]"}]}]}]], "Output", + CellChangeTimes->{ + 3.850657124786903*^9, 3.8506573427936964`*^9, {3.8506573800984526`*^9, + 3.8506574149535217`*^9}, 3.850658116378392*^9, 3.8506581510124197`*^9, + 3.850658191478963*^9}, + CellLabel-> + "Out[135]=",ExpressionUUID->"030cf037-69dd-4ede-945a-154d35f58dfb"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Vertex1", "=", + RowBox[{"{", + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "'"}], "[", "0", "]"}], "\[Equal]", "0"}], + "}"}]}]], "Input", + CellChangeTimes->{{3.834678964818534*^9, 3.834679054548163*^9}, { + 3.83467909814336*^9, 3.8346791212260847`*^9}, {3.834680117471212*^9, + 3.834680131890923*^9}, 3.8506535259593425`*^9, 3.8506542366069365`*^9, + 3.850654566481948*^9, {3.850654609239808*^9, 3.8506546319391756`*^9}, { + 3.850655754862317*^9, 3.850655755260985*^9}, {3.8506561282878056`*^9, + 3.850656128506816*^9}, {3.850656179326333*^9, 3.850656194903241*^9}, { + 3.8506562948960257`*^9, 3.850656298173747*^9}, {3.8506568182908597`*^9, + 3.8506568253680487`*^9}, {3.85065687610748*^9, 3.8506568891480427`*^9}, { + 3.850657251237673*^9, 3.8506572879239845`*^9}}, + CellLabel-> + "In[136]:=",ExpressionUUID->"2298dc58-9998-4336-8637-67bd3311f999"], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "1"], "\[Prime]", + MultilineFunction->None], "[", "0", "]"}], "\[Equal]", "0"}], + "}"}]], "Output", + CellChangeTimes->{ + 3.85065642381654*^9, 3.85065657920348*^9, 3.8506567054719696`*^9, + 3.8506567758191576`*^9, 3.8506568930209246`*^9, 3.8506569497505836`*^9, + 3.850657008890003*^9, {3.85065709545033*^9, 3.8506571249798975`*^9}, + 3.850657342905113*^9, {3.8506573802064686`*^9, 3.8506574150084963`*^9}, + 3.85065811643744*^9, 3.8506581510644197`*^9, 3.850658191488014*^9}, + CellLabel-> + "Out[136]=",ExpressionUUID->"a6be866c-ba05-46da-a41c-ad1aa50b9437"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"Vertex2", "=", + RowBox[{"{", + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "'"}], "[", "0", "]"}], "-", + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "'"}], "[", "\[Pi]", "]"}], "+", + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "'"}], "[", "2", "]"}], "+", + RowBox[{ + SubscriptBox["y", "1"], "[", "2", "]"}]}], "\[Equal]", "1"}], ",", + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "[", "0", "]"}], "==", + RowBox[{ + SubscriptBox["y", "2"], "[", "\[Pi]", "]"}], "==", + RowBox[{ + SubscriptBox["y", "1"], "[", "2", "]"}]}]}], "}"}]}]], "Input", + CellChangeTimes->{{3.834678815307212*^9, 3.834678960354398*^9}, + 3.834679091815206*^9, {3.834680101997093*^9, 3.834680108768464*^9}, { + 3.8506535281132536`*^9, 3.850653528401059*^9}, 3.850654240775215*^9, { + 3.8506546241912746`*^9, 3.850654625032473*^9}, {3.8506555832915897`*^9, + 3.850655587715314*^9}, {3.850656105505641*^9, 3.8506561059986706`*^9}, + 3.850656150647132*^9, {3.8506562839373245`*^9, 3.850656287632543*^9}, { + 3.850656791161252*^9, 3.850656812949095*^9}, {3.850657249495758*^9, + 3.8506572496053743`*^9}, {3.850657299623186*^9, 3.850657336922967*^9}, { + 3.8506580963518257`*^9, 3.850658107182722*^9}, 3.8506581832011404`*^9}, + CellLabel-> + "In[137]:=",ExpressionUUID->"970db54f-1695-48ab-80b5-d0a38cf1f9d2"], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{ + RowBox[{ + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "[", "2", "]"}], "+", + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "1"], "\[Prime]", + MultilineFunction->None], "[", "2", "]"}], "+", + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "2"], "\[Prime]", + MultilineFunction->None], "[", "0", "]"}], "-", + RowBox[{ + SuperscriptBox[ + SubscriptBox["y", "2"], "\[Prime]", + MultilineFunction->None], "[", "\[Pi]", "]"}]}], "\[Equal]", "1"}], + ",", + RowBox[{ + RowBox[{ + SubscriptBox["y", "2"], "[", "0", "]"}], "\[Equal]", + RowBox[{ + SubscriptBox["y", "2"], "[", "\[Pi]", "]"}], "\[Equal]", + RowBox[{ + SubscriptBox["y", "1"], "[", "2", "]"}]}]}], "}"}]], "Output", + CellChangeTimes->{ + 3.8346789636705008`*^9, 3.834680109930167*^9, 3.85065424687862*^9, + 3.8506546578932705`*^9, 3.850655599313405*^9, 3.850655759382452*^9, + 3.8506562018607187`*^9, 3.8506563089355545`*^9, 3.850656423759699*^9, + 3.850656579146865*^9, 3.850656705414361*^9, 3.8506567757631865`*^9, + 3.850656892961679*^9, 3.8506569496581817`*^9, 3.8506570088279743`*^9, { + 3.8506570953848543`*^9, 3.85065712492*^9}, {3.8506573411184406`*^9, + 3.8506573429656816`*^9}, {3.850657380269035*^9, 3.8506574150704403`*^9}, + 3.8506581164970045`*^9, 3.850658151120613*^9, 3.850658191557211*^9}, + CellLabel-> + "Out[137]=",ExpressionUUID->"1298c15d-68d8-4fc9-8a92-0888690e3cad"] +}, Open ]], + +Cell[CellGroupData[{ + +Cell[BoxData[ + RowBox[{"DSolve", "[", + RowBox[{ + RowBox[{"{", + RowBox[{"Side1", ",", "Side2", ",", "Vertex1", ",", "Vertex2"}], "}"}], + ",", + RowBox[{"{", + RowBox[{ + SubscriptBox["y", "1"], ",", + SubscriptBox["y", "2"]}], "}"}], ",", "x"}], "]"}]], "Input", + CellChangeTimes->{{3.834679159495409*^9, 3.834679199735717*^9}, { + 3.85065464029145*^9, 3.850654651576705*^9}}, + CellLabel-> + "In[138]:=",ExpressionUUID->"2199c3ab-f46c-4b40-813c-114a1efe8a45"], + +Cell[BoxData[ + RowBox[{"{", + RowBox[{"{", + RowBox[{ + RowBox[{ + SubscriptBox["y", "1"], "\[Rule]", + RowBox[{"Function", "[", + RowBox[{ + RowBox[{"{", "x", "}"}], ",", + RowBox[{"Cos", "[", + RowBox[{"\[Pi]", " ", "x"}], "]"}]}], "]"}]}], ",", + RowBox[{ + SubscriptBox["y", "2"], "\[Rule]", + RowBox[{"Function", "[", + RowBox[{ + RowBox[{"{", "x", "}"}], ",", + RowBox[{"Cos", "[", + RowBox[{"2", " ", "x"}], "]"}]}], "]"}]}]}], "}"}], "}"}]], "Output", + CellChangeTimes->{ + 3.8346792010932693`*^9, 3.8346801435371943`*^9, 3.8506542471161647`*^9, + 3.8506546580415635`*^9, 3.8506555994399376`*^9, 3.850655759509415*^9, + 3.8506562019967995`*^9, 3.850656309179753*^9, 3.850656423951334*^9, + 3.85065657934756*^9, 3.8506567056216993`*^9, 3.8506567759728107`*^9, + 3.8506568930810633`*^9, 3.8506569498532786`*^9, 3.8506570089621677`*^9, { + 3.850657095533163*^9, 3.850657125050787*^9}, 3.850657343041707*^9, { + 3.8506573803445215`*^9, 3.8506574151474514`*^9}, 3.8506581165773277`*^9, + 3.850658151195318*^9, 3.8506581916392717`*^9}, + CellLabel-> + "Out[138]=",ExpressionUUID->"9c3e8695-8884-4567-80de-a5e222a16a12"] +}, Open ]] +}, Open ]] +}, +WindowSize->{1016.25, 622.5}, +WindowMargins->{{0, Automatic}, {Automatic, 0}}, +FrontEndVersion->"12.1 for Microsoft Windows (64-bit) (June 19, 2020)", +StyleDefinitions->"Default.nb", +ExpressionUUID->"248e3ec3-78ca-428e-ab73-ca0173c7ad94" +] +(* End of Notebook Content *) + +(* Internal cache information *) +(*CellTagsOutline +CellTagsIndex->{} +*) +(*CellTagsIndex +CellTagsIndex->{} +*) +(*NotebookFileOutline +Notebook[{ +Cell[CellGroupData[{ +Cell[580, 22, 314, 8, 54, "Subsection",ExpressionUUID->"459730e3-b4c0-4dc3-86d6-d1c4affa7b98"], +Cell[CellGroupData[{ +Cell[919, 34, 665, 16, 28, "Input",ExpressionUUID->"c2ce4243-f67f-433e-af2d-d14d5a6e465c"], +Cell[1587, 52, 918, 20, 32, "Output",ExpressionUUID->"a0571591-7e57-47a7-8111-76c01f353f9a"] +}, Open ]], +Cell[CellGroupData[{ +Cell[2542, 77, 681, 16, 28, "Input",ExpressionUUID->"b1e599e2-c453-4afb-b14b-eb0509d45b0e"], +Cell[3226, 95, 504, 14, 32, "Output",ExpressionUUID->"b1330e79-85cb-4ddf-9406-88895e4a7d40"] +}, Open ]], +Cell[CellGroupData[{ +Cell[3767, 114, 915, 18, 28, "Input",ExpressionUUID->"2ca5a43f-0e7c-4d8f-a01a-3d97f0f7c774"], +Cell[4685, 134, 623, 15, 32, "Output",ExpressionUUID->"a697f319-b098-49f5-afe1-f7eb7fc72992"] +}, Open ]], +Cell[CellGroupData[{ +Cell[5345, 154, 1309, 32, 28, "Input",ExpressionUUID->"fbbf2e80-0c95-4928-bfff-f47af7a9fe52"], +Cell[6657, 188, 1380, 35, 32, "Output",ExpressionUUID->"782125e5-9d48-4f0c-85a8-170051501035"] +}, Open ]], +Cell[CellGroupData[{ +Cell[8074, 228, 481, 13, 28, "Input",ExpressionUUID->"f6d4efcc-c2fd-430f-a82d-c896e07895ac"], +Cell[8558, 243, 832, 16, 42, "Message",ExpressionUUID->"fba1ea3b-1021-40d6-8754-0dd930352e64"], +Cell[9393, 261, 1287, 34, 32, "Output",ExpressionUUID->"7dc5a94a-4de5-4e7c-be83-fcc6e20b3b3d"] +}, Open ]] +}, Open ]], +Cell[CellGroupData[{ +Cell[10729, 301, 369, 9, 54, "Subsection",ExpressionUUID->"6ffa0d99-169a-43b4-89fe-f3d442d8d5fc"], +Cell[CellGroupData[{ +Cell[11123, 314, 665, 16, 28, "Input",ExpressionUUID->"a25ca004-781d-4237-bb37-d5989eb275ed"], +Cell[11791, 332, 966, 21, 32, "Output",ExpressionUUID->"41e24e3d-6bb6-4743-b3b4-e1702e2bbd7a"] +}, Open ]], +Cell[CellGroupData[{ +Cell[12794, 358, 681, 16, 28, "Input",ExpressionUUID->"1a990e1b-45f5-4a75-b891-c354571430b0"], +Cell[13478, 376, 554, 15, 32, "Output",ExpressionUUID->"030cf037-69dd-4ede-945a-154d35f58dfb"] +}, Open ]], +Cell[CellGroupData[{ +Cell[14069, 396, 915, 18, 28, "Input",ExpressionUUID->"2298dc58-9998-4336-8637-67bd3311f999"], +Cell[14987, 416, 668, 15, 32, "Output",ExpressionUUID->"a6be866c-ba05-46da-a41c-ad1aa50b9437"] +}, Open ]], +Cell[CellGroupData[{ +Cell[15692, 436, 1449, 34, 28, "Input",ExpressionUUID->"970db54f-1695-48ab-80b5-d0a38cf1f9d2"], +Cell[17144, 472, 1493, 37, 32, "Output",ExpressionUUID->"1298c15d-68d8-4fc9-8a92-0888690e3cad"] +}, Open ]], +Cell[CellGroupData[{ +Cell[18674, 514, 481, 13, 28, "Input",ExpressionUUID->"2199c3ab-f46c-4b40-813c-114a1efe8a45"], +Cell[19158, 529, 1200, 28, 81, "Output",ExpressionUUID->"9c3e8695-8884-4567-80de-a5e222a16a12"] +}, Open ]] +}, Open ]] +} +] +*) + +(* End of internal cache information *) + diff --git a/source/examples/poissonExampleLollipopHomogeneous.mlx b/source/examples/poissonExampleLollipopHomogeneous.mlx new file mode 100644 index 00000000..1b2f64b8 Binary files /dev/null and b/source/examples/poissonExampleLollipopHomogeneous.mlx differ diff --git a/source/examples/poissonExamplePlotCoords.m b/source/examples/poissonExamplePlotCoords.m new file mode 100644 index 00000000..78c5d184 --- /dev/null +++ b/source/examples/poissonExamplePlotCoords.m @@ -0,0 +1,28 @@ +function plotCoords = poissonExamplePlotCoords(Phi) + +LVec=Phi.L; + +plotCoords.x1Node=[0;4;4]; +plotCoords.x2Node=[0;0;-1]; + +plotCoords.x1Edge = cell(4,1); +plotCoords.x2Edge = cell(4,1); + +% Edge 1 +theta=Phi.distributePlotCoords(0,2*pi,1); +plotCoords.x1Edge{1} = cos(theta)-1; +plotCoords.x2Edge{1} = sin(theta); + +% Edge 2 +plotCoords.x1Edge{2} = Phi.distributePlotCoords(0,4,2); +plotCoords.x2Edge{2} = zeros(size(plotCoords.x1Edge{2})); + +% Edge 3 +theta=Phi.distributePlotCoords(0,2*pi,3); +plotCoords.x1Edge{3} = sin(theta)+4; +plotCoords.x2Edge{3} = -cos(theta)+1; + +% Edge 4 +plotCoords.x2Edge{4} = Phi.distributePlotCoords(0,-1,4); +plotCoords.x1Edge{4} = 4+zeros(size(plotCoords.x2Edge{4})); + diff --git a/source/examples/starEigenfuctionsDemo.mlx b/source/examples/starEigenfuctionsDemo.mlx index 402a4f4c..9911a89f 100644 Binary files a/source/examples/starEigenfuctionsDemo.mlx and b/source/examples/starEigenfuctionsDemo.mlx differ diff --git a/source/examples/time evolution examples/dumbbellHeatDemo.mlx b/source/examples/time evolution examples/dumbbellHeatDemo.mlx new file mode 100644 index 00000000..8030e818 Binary files /dev/null and b/source/examples/time evolution examples/dumbbellHeatDemo.mlx differ diff --git a/source/examples/time evolution examples/dumbbellNLSDemo.mlx b/source/examples/time evolution examples/dumbbellNLSDemo.mlx new file mode 100644 index 00000000..fed15acb Binary files /dev/null and b/source/examples/time evolution examples/dumbbellNLSDemo.mlx differ diff --git a/source/examples/time evolution examples/dumbbellNLSDemo1.mlx b/source/examples/time evolution examples/dumbbellNLSDemo1.mlx new file mode 100644 index 00000000..02dc59a6 Binary files /dev/null and b/source/examples/time evolution examples/dumbbellNLSDemo1.mlx differ diff --git a/source/examples/time evolution examples/dumbbellNLSDemo2.mlx b/source/examples/time evolution examples/dumbbellNLSDemo2.mlx new file mode 100644 index 00000000..0cb18032 Binary files /dev/null and b/source/examples/time evolution examples/dumbbellNLSDemo2.mlx differ diff --git a/source/examples/time evolution examples/dumbbellNLSTimePeriodic.mlx b/source/examples/time evolution examples/dumbbellNLSTimePeriodic.mlx new file mode 100644 index 00000000..d13e23dc Binary files /dev/null and b/source/examples/time evolution examples/dumbbellNLSTimePeriodic.mlx differ diff --git a/source/examples/time evolution examples/dumbbellNLSUniformDemo.mlx b/source/examples/time evolution examples/dumbbellNLSUniformDemo.mlx new file mode 100644 index 00000000..b84d6683 Binary files /dev/null and b/source/examples/time evolution examples/dumbbellNLSUniformDemo.mlx differ diff --git a/source/examples/time evolution examples/dumbbellTimeEvolution4DAEsDemo.mlx b/source/examples/time evolution examples/dumbbellTimeEvolution4DAEsDemo.mlx deleted file mode 100644 index 765160b9..00000000 Binary files a/source/examples/time evolution examples/dumbbellTimeEvolution4DAEsDemo.mlx and /dev/null differ diff --git a/source/examples/time evolution examples/lineTimeEvolution4DAEsDemo.mlx b/source/examples/time evolution examples/lineDemo.mlx similarity index 100% rename from source/examples/time evolution examples/lineTimeEvolution4DAEsDemo.mlx rename to source/examples/time evolution examples/lineDemo.mlx diff --git a/source/examples/time evolution examples/lollipopSchrodingerDemo.mlx b/source/examples/time evolution examples/lollipopSchrodingerDemo.mlx new file mode 100644 index 00000000..152c9ce6 Binary files /dev/null and b/source/examples/time evolution examples/lollipopSchrodingerDemo.mlx differ diff --git a/source/examples/time evolution examples/starTimeEvolution4DAEsDemo.mlx b/source/examples/time evolution examples/starTimeEvolution4DAEsDemo.mlx deleted file mode 100644 index 260fe17e..00000000 Binary files a/source/examples/time evolution examples/starTimeEvolution4DAEsDemo.mlx and /dev/null differ diff --git a/source/examples/time evolution examples/starTimeEvolutionDemo.mlx b/source/examples/time evolution examples/starTimeEvolutionDemo.mlx new file mode 100644 index 00000000..943933d1 Binary files /dev/null and b/source/examples/time evolution examples/starTimeEvolutionDemo.mlx differ diff --git a/source/templates/BMP.m b/source/templates/BMP.m index ae36b027..28f9dce9 100644 --- a/source/templates/BMP.m +++ b/source/templates/BMP.m @@ -4,6 +4,7 @@ opts.LVec {mustBeNumeric} = [2 pi pi pi 2]; opts.nX {mustBeNumeric} = 10; opts.robinCoeff {mustBeNumeric} = 0; + opts.Discretization {mustBeNonzeroLengthText, mustBeMember(opts.Discretization,{'Uniform','Chebyshev'})} = 'Uniform'; end assert( any(length(opts.nX)==[1 2 5]), 'nX must have length 1, 2, or 5') diff --git a/source/templates/multibell.m b/source/templates/multibell.m index 41f33317..a7b1ad19 100644 --- a/source/templates/multibell.m +++ b/source/templates/multibell.m @@ -5,7 +5,7 @@ % nBells = 1 returns the lollipop % nBells = 2 returns the dumbbell % nBells >= 3 returns the multibell, which consists of lollipops joined at -% the "bottoms of the sticks" +% the "bottoms of the sticks" arguments opts.LVec {mustBeNumeric} = [2 pi]; diff --git a/source/templates/nbell.m b/source/templates/nbell.m index 2ab30aca..7f2f0f9e 100644 --- a/source/templates/nbell.m +++ b/source/templates/nbell.m @@ -6,6 +6,7 @@ opts.nX {mustBeNumeric} = 20; opts.robinCoeff {mustBeNumeric} = 0; opts.Discretization {mustBeNonzeroLengthText, mustBeMember(opts.Discretization,{'Uniform','Chebyshev'})} = 'Uniform'; + opts.nodeData {mustBeNumeric} = [0 0]; end assert(mod(length(opts.LVec),2)==0,'The length vector must have an even number of elements') @@ -33,4 +34,4 @@ target=[nPlus nPlus]; Phi = quantumGraph(source, target,opts.LVec,'nxVec',opts.nX,'robinCoeff',opts.robinCoeff,... - 'Discretization',opts.Discretization); \ No newline at end of file + 'Discretization',opts.Discretization,'nodeData',opts.nodeData); \ No newline at end of file diff --git a/source/utilities/myMenu.m b/source/utilities/myMenu.m deleted file mode 100644 index c8ad66a7..00000000 --- a/source/utilities/myMenu.m +++ /dev/null @@ -1,299 +0,0 @@ -function k = myMenu(xHeader,varargin) -%MENU Generate a menu of choices for user input. -% CHOICE = MENU(HEADER, ITEM1, ITEM2, ... ) displays the HEADER -% string followed in sequence by the menu-item strings: ITEM1, ITEM2, -% ... ITEMn. Returns the number of the selected menu-item as CHOICE, -% a scalar value. There is no limit to the number of menu items. -% -% CHOICE = MENU(HEADER, ITEMLIST) where ITEMLIST is a string, cell -% array is also a valid syntax. -% -% On most graphics terminals MENU will display the menu-items as push -% buttons in a figure window, otherwise they will be given as a numbered -% list in the command window (see example, below). -% -% Example: -% K = menu('Choose a color','Red','Blue','Green') -% %creates a figure with buttons labeled 'Red', 'Blue' and 'Green' -% %The button clicked by the user is returned as K (i.e. K = 2 -% implies that the user selected Blue). -% -% See also UICONTROL, UIMENU, GUIDE. - -% Copyright 1984-2012 The MathWorks, Inc. - -if nargin > 0 - xHeader = convertStringsToChars(xHeader); -end - -if nargin > 1 - [varargin{:}] = convertStringsToChars(varargin{:}); -end - -%========================================================================= -% Check input -%------------------------------------------------------------------------- -if nargin < 2 - disp(getString(message('MATLAB:uistring:menu:NoMenuItemsToChooseFrom'))) - k=0; - return; -elseif nargin==2 && iscell(varargin{1}) - ArgsIn = varargin{1}; % a cell array was passed in -else - ArgsIn = varargin; % use the varargin cell array -end - -%------------------------------------------------------------------------- -% Create the appropriate menu based on HG availability -%------------------------------------------------------------------------- -if matlab.ui.internal.isFigureShowEnabled() - % Create a GUI menu to acquire answer "k" - k = local_GUImenu( xHeader, ArgsIn ); -else - % Create an ascii menu to acquire answer "k" - k = local_ASCIImenu( xHeader, ArgsIn ); -end % if - -%%######################################################################### -% END : main function "menu" -%%######################################################################### - -%%######################################################################### -% BEGIN : local function local_ASCIImenu -%%######################################################################### -function k = local_ASCIImenu( xHeader, xcItems ) - -% local function to display an ascii-generated menu and return the user's -% selection from that menu as an index into the xcItems cell array - -%------------------------------------------------------------------------- -% Calculate the number of items in the menu -%------------------------------------------------------------------------- -numItems = length(xcItems); - -%------------------------------------------------------------------------- -% Continuous loop to redisplay menu until the user makes a valid choice -%------------------------------------------------------------------------- -while 1 - % Display the header - disp(' ') - disp(['----- ',xHeader,' -----']) - disp(' ') - % Display items in a numbered list - for n = 1 : numItems - disp( [ ' ' int2str(n) ') ' xcItems{n} ] ) - end - disp(' ') - % Prompt for user input - k = input('Select a menu number: '); - % Check input: - % 1) make sure k has a value - if isempty(k), k = -1; end - % 2) make sure the value of k is valid - if (k < 1) || (k > numItems) ... - || ~isa(k,'double') ... - || ~isreal(k) || (isnan(k)) || isinf(k) - % Failed a key test. Ask question again - disp(' ') - disp(getString(message('MATLAB:uistring:menu:SelectionOutOfRangeTryAgain'))) - else - % Passed all tests, exit loop and return k - return - end % if k... -end % while 1 - -%%######################################################################### -% END : local function local_ASCIImenu -%%######################################################################### - -%%######################################################################### -% BEGIN : local function local_GUImenu -%%######################################################################### -function k = local_GUImenu( xHeader, xcItems ) - -% local function to display a Handle Graphics menu and return the user's -% selection from that menu as an index into the xcItems cell array - -%========================================================================= -% SET UP -%========================================================================= -% Set spacing and sizing parameters for the GUI elements -%------------------------------------------------------------------------- -MenuUnits = 'pixels'; % units used for all HG objects -textPadding = [22 12]; % extra [Width Height] on uicontrols to pad text -uiGap = 5; % space between uicontrols -uiBorder = 10; % space between edge of figure and any uicontrol -winTopGap = 300; % gap between top of screen and top of figure ** -winLeftGap = 600; % gap between side of screen and side of figure ** -winWideMin = 140; % minimum window width necessary to show title - -% ** "figure" ==> viewable figure. You must allow space for the OS to add -% a title bar (aprx 42 points on Mac and Windows) and a window border -% (usu 2-6 points). Otherwise user cannot move the window. - -%------------------------------------------------------------------------- -% Calculate the number of items in the menu -%------------------------------------------------------------------------- -numItems = length( xcItems ); - -%========================================================================= -% BUILD -%========================================================================= -% Create a generically-sized invisible figure window -%------------------------------------------------------------------------ -menuFig = figure( 'WindowStyle', 'normal', ... - 'Units' ,MenuUnits, ... - 'Visible' ,'off', ... - 'NumberTitle' ,'off', ... - 'Name' ,getString(message('MATLAB:uistring:menu:MENU')), ... - 'Resize' ,'off', ... - 'Colormap' ,[], ... - 'MenuBar' ,'none',... - 'ToolBar' ,'none' ... - ); - -%------------------------------------------------------------------------ -% Add generically-sized header text with same background color as figure -%------------------------------------------------------------------------ -hText = uicontrol(... - 'Parent' ,menuFig, ... - 'Style' ,'text', ... - 'String' ,xHeader, ... - 'Units' ,MenuUnits, ... - 'Position' ,[ 600 500 100 20 ], ... - 'HorizontalAlignment' ,'center',... - 'BackgroundColor' ,get(menuFig,'Color') ); - -% Record extent of text string -maxsize = get( hText, 'Extent' ); -textWide = maxsize(3); -textHigh = maxsize(4); - -%------------------------------------------------------------------------ -% Add generically-spaced buttons below the header text -%------------------------------------------------------------------------ -% Loop to add buttons in reverse order (to automatically initialize numitems). -% Note that buttons may overlap, but are placed in correct position relative -% to each other. They will be resized and spaced evenly later on. - -hBtn = zeros(numItems, 1); -for idx = numItems : -1 : 1 % start from top of screen and go down - n = numItems - idx + 1; % start from 1st button and go to last - % make a button - hBtn(n) = uicontrol( ... - 'Parent' ,menuFig, ... - 'Units' ,MenuUnits, ... - 'Position' ,[uiBorder uiGap*idx textHigh textWide], ... - 'Callback' , {@menucallback, n}, ... - 'String' ,xcItems{n} ); -end % for - -%========================================================================= -% TWEAK -%========================================================================= -% Calculate Optimal UIcontrol dimensions based on max text size -%------------------------------------------------------------------------ -cAllExtents = get( hBtn, {'Extent'} ); % put all data in a cell array -AllExtents = cat( 1, cAllExtents{:} ); % convert to an n x 3 matrix -maxsize = max( AllExtents(:,3:4) ); % calculate the largest width & height -maxsize = maxsize + textPadding; % add some blank space around text -btnHigh = maxsize(2); -btnWide = maxsize(1); - -%------------------------------------------------------------------------ -% Retrieve screen dimensions (in correct units) -%------------------------------------------------------------------------ -screensize = get(0,'ScreenSize'); % record screensize - -%------------------------------------------------------------------------ -% How many rows and columns of buttons will fit in the screen? -% Note: vertical space for buttons is the critical dimension -% --window can't be moved up, but can be moved side-to-side -%------------------------------------------------------------------------ -openSpace = screensize(4) - winTopGap - 2*uiBorder - textHigh; -numRows = min( floor( openSpace/(btnHigh + uiGap) ), numItems ); -if numRows == 0; numRows = 1; end % Trivial case--but very safe to do -numCols = ceil( numItems/numRows ); - -%------------------------------------------------------------------------ -% Resize figure to place it in top left of screen -%------------------------------------------------------------------------ -% Calculate the window size needed to display all buttons -winHigh = numRows*(btnHigh + uiGap) + textHigh + 2*uiBorder; -winWide = numCols*(btnWide) + (numCols - 1)*uiGap + 2*uiBorder; - -% Make sure the text header fits -if winWide < (2*uiBorder + textWide) - winWide = 2*uiBorder + textWide; -end - -% Make sure the dialog name can be shown -if winWide < winWideMin %pixels - winWide = winWideMin; -end - -% Determine final placement coordinates for bottom of figure window -bottom = screensize(4) - (winHigh + winTopGap); - -% Set figure window position -set( menuFig, 'Position', [winLeftGap bottom winWide winHigh] ); - -%------------------------------------------------------------------------ -% Size uicontrols to fit everyone in the window and see all text -%------------------------------------------------------------------------ -% Calculate coordinates of bottom-left corner of all buttons -xPos = ( uiBorder + (0:numCols-1)'*( btnWide + uiGap )*ones(1,numRows) )'; -xPos = xPos(1:numItems); % [ all 1st col; all 2nd col; ...; all nth col ] -yPos = ( uiBorder + (numRows-1:-1:0)'*( btnHigh + uiGap )*ones(1,numCols) ); -yPos = yPos(1:numItems); % [ rows 1:m; rows 1:m; ...; rows 1:m ] - -% Combine with desired button size to get a cell array of position vectors -allBtn = ones(numItems,1); -uiPosMtx = [ xPos(:), yPos(:), btnWide*allBtn, btnHigh*allBtn ]; -cUIPos = num2cell( uiPosMtx( 1:numItems, : ), 2 ); - -% adjust all buttons -set( hBtn, {'Position'}, cUIPos ); - -%------------------------------------------------------------------------ -% Align the Text and Buttons horizontally and distribute them vertically -%------------------------------------------------------------------------ - -% Calculate placement position of the Header -textWide = winWide - 2*uiBorder; - -% Move Header text into correct position near the top of figure -set( hText, ... - 'Position', [ uiBorder winHigh-uiBorder-textHigh textWide textHigh ] ); - -%========================================================================= -% ACTIVATE -%========================================================================= -% Make figure visible -%------------------------------------------------------------------------ -set( menuFig, 'Visible', 'on' ); - -%------------------------------------------------------------------------ -% Wait for choice to be made (i.e UserData must be assigned)... -%------------------------------------------------------------------------ -waitfor(menuFig,'userdata') - -%------------------------------------------------------------------------ -% Selection has been made or figure has been deleted. -% Assign k and delete the Menu figure if it is still valid. -%------------------------------------------------------------------------ -if ishghandle(menuFig) - k = get(menuFig,'UserData'); - delete(menuFig) -else - % The figure was deleted without a selection. Return 0. - k = 0; -end - -%%######################################################################### -% END : local function local_GUImenu -%%######################################################################### - - -function menucallback(btn, evd, index) %#ok -set(gcbf, 'UserData', index);