Fix the standard basis function ek.m

resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolve15s.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/@quantumGraph.type.File/timeEvolveARK.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/continuation.type.File/HHBDetector.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/continuation.type.File/HHBLocator.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info />

resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/1.type.DIR_SIGNIFIER.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info />

resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/HHBDetector.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/HHBLocator.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/adjointContinuationMethodFast.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/adjointContinuationMethodShift.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/continuation.type.File/acm.type.File/eigsJL.m.type.File.xml

@@ -0,0 +1,6 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.Files/source.type.File/examples.type.File/time evolution examples.type.File.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info />

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 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info />

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 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

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 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

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 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

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 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

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 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design" />
+    </Category>
+</Info>

resources/project/Root.type.ProjectPath/d4264395-e014-41c5-8cf6-74b2dd45f9a9.type.Reference.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info Ref="source/continuation/acm" Type="Relative" />

resources/project/Root.type.ProjectPath/e6a4ab36-b927-44c8-9c3b-d8c041cd946a.type.Reference.xml

@@ -0,0 +1,2 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<Info Ref="source/examples/time evolution examples" Type="Relative" />

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

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

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;
 
 

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;

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

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;

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

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

source/@quantumGraph/quantumGraph.m

@@ -6,8 +6,10 @@
     properties
         qg      % the main quantum graph, a digraph
         discretization;
+        derivative;
         laplacianMatrix;
         weightMatrix;
+        weightMatrixWithBCs;
         vertexConditionAssignmentMatrix;
         explicitLaplacian;
     end

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

source/@quantumGraph/secularDet.m

@@ -28,15 +28,15 @@
     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;
                 else
                     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);

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