Skip to content

Commit ffdf081

Browse files
author
SPearce
committed
Support for interface problems where the two systems are different sizes on each side
1 parent 07594f1 commit ffdf081

File tree

3 files changed

+35
-41
lines changed

3 files changed

+35
-41
lines changed

‎CompoundMatrixMethod-0.91.paclet‎

10.1 KB
Binary file not shown.

‎CompoundMatrixMethod/CompoundMatrixMethod.m‎

Lines changed: 34 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -5,11 +5,11 @@
55
(* :Title: CompoundMatrixMethod *)
66
(* :Author: Simon Pearce <[email protected]> *)
77
(* :Context: CompoundMatrixMethod` *)
8-
(* :Version: 0.9 *)
9-
(* :Date: 2018-09-26 *)
8+
(* :Version: 0.91 *)
9+
(* :Date: 2022-02-09 *)
1010

1111
(* :Mathematica Version: 10+ *)
12-
(* :Copyright: (c) 2017-18 Simon Pearce *)
12+
(* :Copyright: (c) 2017-22 Simon Pearce *)
1313

1414
BeginPackage["CompoundMatrixMethod`"];
1515

@@ -137,20 +137,20 @@ with radius controlled by ContourRadius (default 1), for a system defined from T
137137
Transpose[Table[Coefficient[minorsTab, \[FormalPhi][i]], {i, 1, Binomial[len, len2]}]]
138138
]
139139

140-
rulesFG[len_?NumericQ] := rulesFG[len] = {\[FormalPhi]\[FormalPhi]L[{l_}][x_] :> Sum[F[l, a] \[FormalPhi]L[a][x], {a, 1, len}],
141-
\[FormalPhi]\[FormalPhi]L[{l_, m_}][x_] :> Sum[F[l, a] F[m, b] \[FormalPhi]L[{a, b}][x], {a, 1, len}, {b, 1, len}],
142-
\[FormalPhi]\[FormalPhi]L[{l_, m_, n_}][x_] :> Sum[F[l, a] F[m, b] F[n, c] \[FormalPhi]L[{a, b, c}][x], {a, 1, len}, {b, 1, len}, {c, 1, len}],
143-
\[FormalPhi]\[FormalPhi]L[{l_, m_, n_, o_}][x_] :> Sum[F[l, a] F[m, b] F[n, c] F[o, d]
144-
\[FormalPhi]L[{a, b, c, d}][x], {a, 1, len}, {b, 1, len}, {c, 1, len}, {d, 1, len}],
145-
\[FormalPhi]\[FormalPhi]L[{l_, m_, n_, o_, p_}][x_] :> Sum[F[l, a] F[m, b] F[n, c] F[o, d] F[p,e]
146-
\[FormalPhi]L[{a, b, c, d, e}][x], {a, 1, len}, {b, 1, len}, {c, 1, len}, {d, 1, len}, {e, 1, len}],
147-
\[FormalPhi]\[FormalPhi]R[{l_}][x_] :> Sum[G[l, a] \[FormalPhi]R[a][x], {a, 1, len}],
148-
\[FormalPhi]\[FormalPhi]R[{l_, m_}][x_] :> Sum[G[l, a] G[m, b] \[FormalPhi]R[{a, b}][x], {a, 1, len}, {b, 1, len}],
149-
\[FormalPhi]\[FormalPhi]R[{l_, m_, n_}][x_] :> Sum[G[l, a] G[m, b] G[n, c] \[FormalPhi]R[{a, b, c}][x], {a, 1, len}, {b, 1, len}, {c, 1, len}],
150-
\[FormalPhi]\[FormalPhi]R[{l_, m_, n_, o_}][x_] :> Sum[G[l, a] G[m, b] G[n, c] G[o, d] \[FormalPhi]R[{a, b, c, d}][
151-
x], {a, 1, len}, {b, 1, len}, {c, 1, len}, {d, 1, len}],
152-
\[FormalPhi]\[FormalPhi]R[{l_, m_, n_, o_, p_}][x_] :> Sum[G[l, a] G[m, b] G[n, c] G[o, d] G[p, e]
153-
\[FormalPhi]R[{a, b, c, d, e}][x], {a, 1, len}, {b, 1, len}, {c, 1, len}, {d, 1, len}, {e, 1, len}]};
140+
rulesF[len_?NumericQ]:= rulesF[len]= {
141+
\[FormalPhi]\[FormalPhi]L[{l_}][x_]:>Sum[F[l,a] \[FormalPhi]L[a][x],{a,1,len}],
142+
\[FormalPhi]\[FormalPhi]L[{l_,m_}][x_]:>Sum[F[l,a] F[m,b] \[FormalPhi]L[{a,b}][x],{a,1,len},{b,1,len}],
143+
\[FormalPhi]\[FormalPhi]L[{l_,m_,n_}][x_]:>Sum[F[l,a] F[m,b] F[n,c] \[FormalPhi]L[{a,b,c}][x],{a,1,len},{b,1,len},{c,1,len}],
144+
\[FormalPhi]\[FormalPhi]L[{l_,m_,n_,o_}][x_]:>Sum[F[l,a] F[m,b] F[n,c] F[o,d] \[FormalPhi]L[{a,b,c,d}][x],{a,1,len},{b,1,len},{c,1,len},{d,1,len}],
145+
\[FormalPhi]\[FormalPhi]L[{l_,m_,n_,o_,p_}][x_]:>Sum[F[l,a] F[m,b] F[n,c] F[o,d] F[p,e] \[FormalPhi]L[{a,b,c,d,e}][x],{a,1,len},{b,1,len},{c,1,len},{d,1,len},{e,1,len}]};
146+
147+
rulesG[len_?NumericQ]:=rulesG[len]={
148+
\[FormalPhi]\[FormalPhi]R[{l_}][x_]:>Sum[G[l,a] \[FormalPhi]R[a][x],{a,1,len}],
149+
\[FormalPhi]\[FormalPhi]R[{l_,m_}][x_]:>Sum[G[l,a] G[m,b] \[FormalPhi]R[{a,b}][x],{a,1,len},{b,1,len}],
150+
\[FormalPhi]\[FormalPhi]R[{l_,m_,n_}][x_]:>Sum[G[l,a] G[m,b] G[n,c] \[FormalPhi]R[{a,b,c}][x],{a,1,len},{b,1,len},{c,1,len}],
151+
\[FormalPhi]\[FormalPhi]R[{l_,m_,n_,o_}][x_]:>Sum[G[l,a] G[m,b] G[n,c] G[o,d] \[FormalPhi]R[{a,b,c,d}][x],{a,1,len},{b,1,len},{c,1,len},{d,1,len}],
152+
\[FormalPhi]\[FormalPhi]R[{l_,m_,n_,o_,p_}][x_]:>Sum[G[l,a] G[m,b] G[n,c] G[o,d] G[p,e] \[FormalPhi]R[{a,b,c,d,e}][x],{a,1,len},{b,1,len},{c,1,len},{d,1,len},{e,1,len}]};
153+
154154

155155

156156
Options[Evans]={NormalizationConstants -> 1, MaxStepFraction->0.01};
@@ -445,7 +445,7 @@ with radius controlled by ContourRadius (default 1), for a system defined from T
445445
Evans[\[FormalLambda]0_?
446446
NumericQ, {ALeftMatrix_?MatrixQ, ARightMatrix_?MatrixQ}, leftBCMatrix_?MatrixQ, rightBCMatrix_?MatrixQ,
447447
{FMatrix_?MatrixQ, GMatrix_?MatrixQ}, {x_ /; !NumericQ[x], xaa_, xmm_, xbb_},opts:OptionsPattern[{Evans,NDSolve}]] :=
448-
Module[{dettt, len, subsets, newYs, leftYICs, rightYICs, phiLeftVector, phiRightVector, LeftBCSolution, RightBCSolution, yLeft, yRight,
448+
Module[{dettt, lenAL,lenAR, subsets, newYsL,newYsR, leftYICs, rightYICs, phiLeftVector, phiRightVector, LeftBCSolution, RightBCSolution, yLeft, yRight,
449449
phiLeft, phiRight, LeftPositiveEigenvalues, RightNegativeEigenvalues, phiLeftICs, phiRightICs, QQ, solutionFromRight,
450450
solutionFromLeft, det, matchPoint,lenLeft,lenRight,subsetsLeft,subsetsRight,QLeft,QRight,xa,xb,xm},
451451

@@ -478,39 +478,34 @@ with radius controlled by ContourRadius (default 1), for a system defined from T
478478
If[! MatrixQ[ARightMatrix /. x -> xb /. \[FormalLambda] -> \[FormalLambda]0 ,NumericQ],
479479
Message[Evans::nonNumericalMatrix, ARightMatrix, xb]; Return[$Failed]];
480480

481-
482-
If[Length[ARightMatrix] != Length[ALeftMatrix],
483-
Message[Evans::MatrixSizesDiffer, ALeftMatrix, ARightMatrix];Return[$Failed]];
484-
485-
len = Length[ARightMatrix];
486-
487-
If[len>10, Message[Evans::InterfaceTooBig, len];Return[$Failed]];
488-
489-
490-
newYs = Through[Array[\[FormalY], {len}][x]];
481+
lenAL=Length[ALeftMatrix];
482+
lenAR=Length[ARightMatrix];
483+
If[Max[lenAL,lenAR]>10,Message[Evans::InterfaceTooBig,Max[lenAL,lenAR]];Return[$Failed]];
484+
newYsL=Through[Array[\[FormalY],{lenAL}][x]];
485+
newYsR=Through[Array[\[FormalY],{lenAR}][x]];
491486

492487
(*Initial conditions for shooting from the LHS*)
493488
LeftBCSolution =
494-
Quiet@Solve[leftBCMatrix.newYs == 0, newYs];
489+
Quiet@Solve[leftBCMatrix.newYsL == 0, newYsL];
495490
leftYICs =
496491
NullSpace[leftBCMatrix /. x -> xa /. \[FormalLambda] -> \[FormalLambda]0, Method -> "DivisionFreeRowReduction"];
497492
lenLeft = Length[leftYICs];
498-
subsetsLeft = Subsets[Range[len], {lenLeft}];
493+
subsetsLeft = Subsets[Range[lenAL], {lenLeft}];
499494
(*Initial conditions for shooting from the RHS*)
500495
RightBCSolution =
501-
Quiet@Solve[rightBCMatrix.newYs == 0, newYs];
496+
Quiet@Solve[rightBCMatrix.newYsR == 0, newYsR];
502497
rightYICs =
503498
NullSpace[rightBCMatrix /. x -> xb /. \[FormalLambda] -> \[FormalLambda]0, Method -> "DivisionFreeRowReduction"];
504499
lenRight = Length[rightYICs];
505-
subsetsRight = Subsets[Range[len], {lenRight}];
500+
subsetsRight = Subsets[Range[lenAR], {lenRight}];
506501
(*Check the initial conditions for each side are enough*)
507502
If[Length[LeftBCSolution] != 1,
508503
Message[Evans::boundarySolutionFailed, xa];
509504
Return[$Failed]];
510505
If[Length[RightBCSolution] != 1,
511506
Message[Evans::boundarySolutionFailed, xb];
512507
Return[$Failed]];
513-
If[Length[leftYICs] + Length[rightYICs] != len,
508+
If[Length[leftYICs]+Length[rightYICs]!=((lenAL+lenAR)/2),
514509
Message[Evans::boundaryConditionRank];
515510
Return[$Failed]];
516511
(*Generate two sets of Phi vaiables,these will be the matrix minors*)
@@ -521,9 +516,9 @@ with radius controlled by ContourRadius (default 1), for a system defined from T
521516
Table[\[FormalPhi]R[i][x], {i, 1, Length[subsetsRight]}];
522517
(*Full set of Initial Conditions for the left and right sides, with the BCs incorporated*)
523518
yLeft =
524-
Transpose[leftYICs + Table[newYs, {lenLeft}] /. LeftBCSolution[[1]] /.Thread[newYs -> 0]];
519+
Transpose[leftYICs + Table[newYsL, {lenLeft}] /. LeftBCSolution[[1]] /.Thread[newYsL -> 0]];
525520
yRight =
526-
Transpose[rightYICs + Table[newYs, {lenRight}] /. RightBCSolution[[1]] /.Thread[newYs -> 0]];
521+
Transpose[rightYICs + Table[newYsR, {lenRight}] /. RightBCSolution[[1]] /.Thread[newYsR -> 0]];
527522
(*Use the initial conditions on the Y vectors to generate initial conditions for the minors phi*)
528523
phiLeft = (Det[(yLeft /. x -> xa /. \[FormalLambda] -> \[FormalLambda]0)[[#]]] & /@subsetsLeft);
529524
phiRight = (Det[(yRight /. x -> xb /. \[FormalLambda] -> \[FormalLambda]0)[[#]]] & /@subsetsRight);
@@ -545,10 +540,10 @@ with radius controlled by ContourRadius (default 1), for a system defined from T
545540
phiRight];
546541
(*Calculate the Q matrix (phi'=Q phi) for each side*)
547542
QLeft =
548-
qMatrix[len, lenLeft] /. \[FormalCapitalA][i_, j_] :>
543+
qMatrix[lenAL, lenLeft] /. \[FormalCapitalA][i_, j_] :>
549544
ALeftMatrix[[i, j]] /. \[FormalLambda] -> \[FormalLambda]0;
550545
QRight =
551-
qMatrix[len, lenRight] /. \[FormalCapitalA][i_, j_] :>
546+
qMatrix[lenAR, lenRight] /. \[FormalCapitalA][i_, j_] :>
552547
ARightMatrix[[i, j]] /. \[FormalLambda] -> \[FormalLambda]0;
553548

554549
(*Solve for integrating from the left and right*)
@@ -564,11 +559,10 @@ with radius controlled by ContourRadius (default 1), for a system defined from T
564559
(* Now we need to account for the jump conditions at the interface, so instead of the normal determinant it needs
565560
modifying by multiplication by the matrices F and G. *)
566561

567-
det = Total@Table[\[FormalPhi]\[FormalPhi]L[i][x]
568-
\[FormalPhi]\[FormalPhi]R[Complement[Range[len], i]][x] (-1)^(Total[Range[lenLeft] + i]), {i, subsetsLeft}];
562+
det = Total@Table[\[FormalPhi]\[FormalPhi]L[i][x] \[FormalPhi]\[FormalPhi]R[Complement[Range[Length@FMatrix],i]][x] (-1)^(Total[Range[lenLeft]+i]),{i,Subsets[Range[Length@FMatrix],{lenLeft}]}];
569563

570564
dettt =
571-
det /.rulesFG[len]/.reprules2 /.Thread[subsetsLeft -> Range[Length[subsetsLeft]]]
565+
det /.rulesF[lenAL]/.rulesG[lenAR]/.reprules2 /.Thread[subsetsLeft -> Range[Length[subsetsLeft]]]
572566
/.Thread[subsetsRight -> Range[Length[subsetsRight]]];
573567

574568
Exp[-Integrate[

‎CompoundMatrixMethod/PacletInfo.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
(* ::Input::Initialization:: *)
44
Paclet[
55
Name->"CompoundMatrixMethod",
6-
Version->"0.9",
6+
Version->"0.91",
77
MathematicaVersion->"10+",
88
Description->"Solve Eigenvalue Boundary Value Problems using the Compound Matrix Method to generate the Evans function. ",
99
Creator->"Simon Pearce",

0 commit comments

Comments
 (0)