9
$\begingroup$

I have a 10x10 linearized BVP which I can write as $$\mathbf{y}'(x) = \mathbf{A}(\omega) \mathbf{y}(x)$$ subject to boundary conditions $$\mathbf{B} \cdot \mathbf{y} = \mathbf{0}, \quad x=0 \\ \mathbf{C} \cdot \mathbf{y} = \mathbf{0}, \quad x=1,$$ where $\mathbf{B}, \mathbf{C}$ are diagonal matrices.

Here is one such value of $A$ (with values rounded), although this depends on a set of other parameters that I'd like to explore how the eigenvalues depend on (particularly looking for stability)

A = SparseArray[{{{1, 2}, {1, 7}, {2, 3}, {3, 4}, {4, 3}, {4, 6}, {5, 6}, {6, 5}, {7, 8}, {9, 10}} -> 1, 
    {{1, 9}, {6, 4}} -> -1, {4, 1} -> -600 ω, {{4, 2}, {4, 7}, {8, 2}} -> -300, 
    {{4, 9}, {8,9}, {10, 2},{10, 7}} -> 300, {8, 7} ->  10000 + 700 ω, {{8, 8}, {10, 10}} -> 14, {10, 9} -> -700 + 700 ω}]

yvec = Through[Array[y, 10][t]];
bcmat = DiagonalMatrix[{0, 1, 1, 0, 1, 0, 0, 1, 0, 1}];
bcs = Thread[Select[Join[bcmat.yvec /.t->0, bcmat.yvec /.t->1],LeafCount[#] > 1 &] == 0]

I'm looking for the (possibly complex) eigenvalues of $\omega$ to give a non-trivial solution to my boundary-value problem, but I get errors regarding numerical ill-conditioning when just trying to use NDSolve, e.g.

NDSolve[Flatten@{Thread[D[yvec, t] == A.yvec]/.ω->1, bcs}, Array[y, 10], {t, 0, 1}]
Improve this question
$\endgroup$
1
  • $\begingroup$ Also, is there no tag for BVPs? $\endgroup$
    – SPPearce
    Sep 5, 2017 at 9:38

1 Answer 1

Reset to default
18
$\begingroup$

Update: This implementation is now a package called CompoundMatrixMethod, hosted on github. It can be installed easily by evaluating:

Needs["PacletManager`"] 
PacletInstall["CompoundMatrixMethod", "Site" -> 
 "http://raw.githubusercontent.com/paclets/Repository/master"]

This version also includes a function ToMatrixSystem which converts a system of ODEs to matrix form (and linearises if necessary), including the boundary conditions. This eliminates the need to set the matrices directly, and also specifies which variable is the eigenvalue, simplifying the notation. Please use the package rather than the code below.

I've written an implementation of the Compound Matrix Method that suits my purposes, and so I'll put it here for other people. A good explanation of this method are available here. Basically the Compound Matrix Method takes an $n$ by $n$ eigenvalue problem of the form $$\mathbf{y}' = A(x, \lambda) \mathbf{y}, \quad a \leq x \leq b, \\ B(x,\lambda) \mathbf{y} = \mathbf{0}, \quad x=a, \\ C(x,\lambda) \mathbf{y} = \mathbf{0}, \quad x=b,$$ and converts it to a larger system of determinants that satisfy a different matrix equation $$ \mathbf{\phi}' = Q(x, \lambda) \mathbf{\phi}.$$ This removes a lot of the stiffness from the equations, as well as being able to also remove the exponential growth terms that dominate away from an eigenvalue.

The code is written for general size $n$, and I've used it for $n=10$. The first time you run the code for a particular size $n$ the general form of matrix $\mathbf{Q}$ will be calculated, for $n=10$ this takes about 3 minutes for me, after that the matrix will be cached. The matching should be independent of the choice of matching point, but you can change it in the code to check that.

reprules = ϕ[a_List] :> Signature[a] ϕ[Sort[a]];
minorsDerivs[list_?VectorQ,len_?NumericQ] := 
 Sum[Sum[AA[y, z] ϕ[list /. y -> z], {z, Union[Complement[Range[len], list], {y}]}], {y, list}] /. reprules

qComponents[n_?NumericQ, len_?NumericQ] := qComponents[n, len] = 
       Coefficient[Table[minorsDerivs[ii, len], {ii, Subsets[Range[len], {len/2}]}] 
       /.Thread[Subsets[Range[len], {len/2}] -> Range[Binomial[len, len/2]]], \[Phi][n]]

Evans[{λ_/;!NumericQ[λ], λλ_?NumericQ}, Amat_?MatrixQ, bvec_?MatrixQ, cvec_?MatrixQ, 
        {x_ /;!NumericQ[x], xa_?NumericQ, xb_?NumericQ,xmatch_:False}] := 
 Module[{ya, yb, ϕpa, ϕmb, valsleft, valsright, ϕpainit, ϕmbinit, posint, 
   negint, ϕmvec, ϕpvec, det, QQ, len, subsets,matchpt},
  len = Length[Amat];

  If[(xa <= xmatch <= xb && NumericQ[xmatch]), matchpt = xmatch, matchpt = (xb - xa)/2];
  If[!EvenQ[len], Print["Matrix A does not have even dimension"]; Abort[]];
  If[Length[Amat] != Length[Transpose[Amat]],Print["Matrix A is not a square matrix"]; Abort[]];
  subsets = Subsets[Range[len], {len/2}];
  ya = NullSpace[bvec];
  If[Length[ya] != len/2, Print["Rank of matrix B is not correct"];Abort[]];
  yb = NullSpace[cvec];
  If[Length[yb] != len/2, Print["Rank of matrix C is not correct"];Abort[]];
  ϕmvec = Table[ϕm[i][x], {i, 1, Length[subsets]}]; 
  ϕpvec = Table[ϕp[i][x], {i, 1, Length[subsets]}];
  ϕpa = (Det[Transpose[ya][[#]]] & /@ subsets);
  ϕmb = (Det[Transpose[yb][[#]]] & /@ subsets);
  valsleft =  Select[Eigenvalues[Amat /. x -> xa /. λ -> λλ], Re[#] > 0 &];
  valsright = Select[Eigenvalues[Amat /. x -> xb /. λ -> λλ], Re[#] < 0 &];
  ϕpainit = Thread[Through[Array[ϕp, {Length[subsets]}][xa]] == ϕpa];
  ϕmbinit = Thread[Through[Array[ϕm, {Length[subsets]}][xb]] == ϕmb];
  QQ = Transpose[Table[qComponents[i, len], {i, 1, Length[subsets]}]] /. 
     AA[i_, j_] :> Amat[[i, j]] /. λ -> λλ;
  posint = NDSolve[{Thread[D[ϕpvec,x] == (QQ - Total[Re@valsleft] IdentityMatrix[Length[QQ]]).ϕpvec], ϕpainit}, 
     Array[ϕp, {Length[subsets]}], {x, xa, xb}][[1]];
  negint = NDSolve[{Thread[D[ϕmvec,x] == (QQ - Total[Re@valsright] IdentityMatrix[Length[QQ]]).ϕmvec], ϕmbinit}, 
     Array[ϕm, {Length[subsets]}], {x, xa, xb}][[1]];
  det = Total@Table[ϕm[i][x] ϕp[Complement[Range[len], i]][x] (-1)^(Total[Range[len/2] + i]) //. reprules /. 
      Thread[subsets -> Range[Length[subsets]]], {i, subsets}];
  Exp[-Integrate[Tr[Amat], {x, xa, matchpt}]] det /. x -> matchpt /. posint /. negint]

For a simple 2nd order eigenvalue problem, $y''(x) + \lambda y(x) = 0, y(0)=y(L)=0$, the roots can be found analytically as $n \pi/L, n \in \mathbb{Z}$. Here the matrix $A$ is {{0,1}, {-\[Lambda]^2, 0}}, and the BCs are DiagonalMatrix[{1, 0}]:

Plot[Evans[{λ, λλ}, {{0, 1}, {-λ^2, 0}}, 
     DiagonalMatrix[{1, 0}], DiagonalMatrix[{1, 0}], {x, 0, 2}], {λλ, 0.1, 20}]

Compound Matrix Method example for simple 2nd order eigenvalue problem

Changing the boundary conditions is straight forward, so for a Robin BCs like $y(0)+2y'(0)=0$ the corresponding matrix $B$ would be {{1, 2}, {0, 0}}.

For the first 4th order example in the linked notes $$\epsilon^4 y''''(x) + 2 \epsilon^2 \lambda \frac{d}{dx}\left[\sin(x) \frac{dy}{dx}\right]+y =0, \\ y(0) = y''(0) = y'(\pi/2) = y'''(\pi/2) = 0,$$ the matrices are given by:

A1={{0,1,0,0}, {0,0,1,0}, {0,0,0,1}, {-1/ϵ^4, -2 ω Cos[x]/ϵ^2, -2 ω Sin[x]/ϵ^2, 0}}; 
B1 = DiagonalMatrix[{1,0,1,0}]; C1 = DiagonalMatrix[{0,1,0,1}];

Evans[{ω, 1}, A1 /. ϵ-> 0.1, B1, C1, {x, 0, Pi/2}]
(* -0.650472 *)

And we can then vary the value of $\omega$ to see the roots:

 Plot[Evans[{ω, ωω}, A1 /.ϵ->0.1, B1, C1, {x, 0, Pi/2}], {ωω, 1, 3}]

Plot of the compound matrix method determinant for a 4x4 case

For a 10x10 example similar to my original question (that has positive eigenvalues):

A2 = {{0, 1, 0, 0, 0, 0, 5, 0, -5, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 
0}, {0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {-625 ω, -(125/2), 2, 
0, 0, 3, -300, 0, 300, 0}, {0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 
0, -1.5, 1/2, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, 0, 
0}, {0, -169, 0, 0, 0, 0, 9175 + 694 ω, 0, 811, 0}, {0, 0, 
0, 0, 0, 0, 0, 0, 0, 1}, {0, 672, 0, 0, 0, 0, 3222, 
0, -709 + 694 ω, 0}};
B2 = C2 = DiagonalMatrix[{0, 1, 1, 0, 1, 0, 0, 1, 0, 1}];

Evans[{ω, 1}, A2, B2, C2, {x, 0, 1}]
(* 0.672945 *)

We can plot and see some positive eigenvalues:

 ListPlot[Table[{ωω,Evans[{ω, ωω}, A2, B2, C2, {x, 0, 1}]},{ωω,0.1,1,0.01}]

Plot of the compound matrix method determinant for a 10x10 case

And then FindRoot will find one:

 FindRoot[Evans[{ω, ωω}, A2, B2, C2, {x, 0, 1}],{ωω,0.5}]

The eigenfunctions can be extracted from this method if required, but I haven't coded that here. The subtraction of the dominant growing eigenvalues from $Q$ may not be suitable for all problems, but is really useful when it works. It will also use exact numbers if you give them in the original matrices, so it'll be faster if you give an approximate number.

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ If anyone has comments about for improving my code feel free to let me know. I was thinking that the root specification should probably be moved to the last option rather than the first (like in functions like FindRoot). $\endgroup$
    – SPPearce
    Sep 19, 2017 at 4:59
  • $\begingroup$ I have read the explanation of this method (in the provided link), where the 1st method (determinantal method) is actually a shooting method? $\endgroup$
    – Nobody
    Feb 9, 2020 at 3:56
  • 1
    $\begingroup$ @Nobody, yes, the determinant method is a shooting method. This is a shooting method too, but the compound matrix method lifts the problem from being on a relatively unstable space of dimension $n$ (which is prone to blowing up to infinity), to one where there is more stable. The Allen and Bridges (2002) paper, doi.org/10.1007/s002110100365, is relatively readable regarding that. The other benefit compared to the determinant method is that it subtracts the exponentially growing parts, which means that the value of the Evans function is order 1, rather than say 10^30. $\endgroup$
    – SPPearce
    Feb 9, 2020 at 6:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.