using NDEigensystem to solve the Mathieu equation Announcing the arrival of Valued Associate...

Does a Draconic Bloodline sorcerer's doubled proficiency bonus for Charisma checks against dragons apply to all dragon types or only the chosen one?

Married in secret, can marital status in passport be changed at a later date?

TV series episode where humans nuke aliens before decrypting their message that states they come in peace

Why isn't everyone flabbergasted about Bran's "gift"?

Are these square matrices always diagonalisable?

Could a cockatrice have parasitic embryos?

Marquee sign letters

How would you suggest I follow up with coworkers about our deadline that's today?

What's the difference between using dependency injection with a container and using a service locator?

Is a self contained air-bullet cartridge feasible?

Why doesn't the university give past final exams' answers?

/bin/ls sorts differently than just ls

Was Objective-C really a hindrance to Apple software development?

Why do people think Winterfell crypts is the safest place for women, children and old people?

What is a 'Key' in computer science?

My admission is revoked after accepting the admission offer

What's called a person who works as someone who puts products on shelves in stores?

Philosophers who were composers?

What's parked in Mil Moscow helicopter plant?

Did war bonds have better investment alternatives during WWII?

using NDEigensystem to solve the Mathieu equation

What to do with someone that cheated their way though university and a PhD program?

Processing ADC conversion result: DMA vs Processor Registers

Simulate round-robin tournament draw



using NDEigensystem to solve the Mathieu equation



Announcing the arrival of Valued Associate #679: Cesar Manara
Unicorn Meta Zoo #1: Why another podcast?How to correctly use DSolve when the force is an impulse (dirac delta) and initial conditions are not zeroIndexing of Large Autonomous System of Equations for Use in NDSolveFEM Solution desired for “Plate with orifice” deflection: Application of Boundary Conditions and use of RegionsSolving an ODE using shooting methodHow to rescale the independent variable?Trouble with shooting method for a 4th-order stiff ODEFinding eigenvalues for Laplacian operator for 3D shape with Neumann boundary conditionsHow do you find the eigenvalues of a PDE (Dynamic Euler-Bernoulli beam)?How to evaluate the PDE solution dependent on the `RegionMarkers"?Using NDEigensystem to solve coupled eigenvalue problem












4












$begingroup$


To be able to apply the differentialequation capabilities of Mathematica to my graduate thesis, I am trying to apply NDEigensystem to an eigenproblem whose solution I know, but I am having some trouble doing so.



As a test problem, I am using an algebraic version of the Mathieu equation,



$$(1-zeta^{2})w^{primeprime}-zeta w^{prime}+left(a+2q-4qzeta^{2}right)w=0$$



For this example I set $q=4/3$ and take only the first three eigenpairs:



m = 3; q = 4/3;
op = -(1 - ζ^2) u''[ζ] + ζ u'[ζ] + 2 q (2 ζ^2 - 1) u[ζ];
bc = DirichletCondition[u[ζ] == 0, True];
{λ, fl} = NDEigensystem[{op, bc}, u, {ζ, 0, 1}, m];


I chose the Mathieu equation as a nontrivial example as Mathematica already has a function for it's evaluation:



λt = Table[MathieuCharacteristicB[2 k, q], {k, m}];
flt = Table[With[{j = j},
MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];


The problem is, I do not get the expected eigenvalues!



λ
(* {4.0708, 17.3259, 39.1877} *)
N[λt]
(* {3.85298, 16.0581, 36.0254} *)


And of course, plotting shows that the eigenequation is not satisfied at all:



With[{u = fl[[1]], b = λ[[1]]},
Plot[(1 - ζ^2) u''[ζ] - ζ u'[ζ] + (b + 2 q - 4 q ζ^2) u[ζ], {ζ, 0, 1}]]
With[{u = flt[[1]], b = λt[[1]]},
Plot[(1 - ζ^2) u''[ζ] - ζ u'[ζ] + (b + 2 q - 4 q ζ^2) u[ζ], {ζ, 0, 1}]]


What was wrong with my attempt? If I can get this example to work, I should be able to apply it to my actual, more complicated problem, so any Good Ideas would be welcome.










share|improve this question









New contributor




宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.







$endgroup$

















    4












    $begingroup$


    To be able to apply the differentialequation capabilities of Mathematica to my graduate thesis, I am trying to apply NDEigensystem to an eigenproblem whose solution I know, but I am having some trouble doing so.



    As a test problem, I am using an algebraic version of the Mathieu equation,



    $$(1-zeta^{2})w^{primeprime}-zeta w^{prime}+left(a+2q-4qzeta^{2}right)w=0$$



    For this example I set $q=4/3$ and take only the first three eigenpairs:



    m = 3; q = 4/3;
    op = -(1 - ζ^2) u''[ζ] + ζ u'[ζ] + 2 q (2 ζ^2 - 1) u[ζ];
    bc = DirichletCondition[u[ζ] == 0, True];
    {λ, fl} = NDEigensystem[{op, bc}, u, {ζ, 0, 1}, m];


    I chose the Mathieu equation as a nontrivial example as Mathematica already has a function for it's evaluation:



    λt = Table[MathieuCharacteristicB[2 k, q], {k, m}];
    flt = Table[With[{j = j},
    MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];


    The problem is, I do not get the expected eigenvalues!



    λ
    (* {4.0708, 17.3259, 39.1877} *)
    N[λt]
    (* {3.85298, 16.0581, 36.0254} *)


    And of course, plotting shows that the eigenequation is not satisfied at all:



    With[{u = fl[[1]], b = λ[[1]]},
    Plot[(1 - ζ^2) u''[ζ] - ζ u'[ζ] + (b + 2 q - 4 q ζ^2) u[ζ], {ζ, 0, 1}]]
    With[{u = flt[[1]], b = λt[[1]]},
    Plot[(1 - ζ^2) u''[ζ] - ζ u'[ζ] + (b + 2 q - 4 q ζ^2) u[ζ], {ζ, 0, 1}]]


    What was wrong with my attempt? If I can get this example to work, I should be able to apply it to my actual, more complicated problem, so any Good Ideas would be welcome.










    share|improve this question









    New contributor




    宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
    Check out our Code of Conduct.







    $endgroup$















      4












      4








      4





      $begingroup$


      To be able to apply the differentialequation capabilities of Mathematica to my graduate thesis, I am trying to apply NDEigensystem to an eigenproblem whose solution I know, but I am having some trouble doing so.



      As a test problem, I am using an algebraic version of the Mathieu equation,



      $$(1-zeta^{2})w^{primeprime}-zeta w^{prime}+left(a+2q-4qzeta^{2}right)w=0$$



      For this example I set $q=4/3$ and take only the first three eigenpairs:



      m = 3; q = 4/3;
      op = -(1 - ζ^2) u''[ζ] + ζ u'[ζ] + 2 q (2 ζ^2 - 1) u[ζ];
      bc = DirichletCondition[u[ζ] == 0, True];
      {λ, fl} = NDEigensystem[{op, bc}, u, {ζ, 0, 1}, m];


      I chose the Mathieu equation as a nontrivial example as Mathematica already has a function for it's evaluation:



      λt = Table[MathieuCharacteristicB[2 k, q], {k, m}];
      flt = Table[With[{j = j},
      MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];


      The problem is, I do not get the expected eigenvalues!



      λ
      (* {4.0708, 17.3259, 39.1877} *)
      N[λt]
      (* {3.85298, 16.0581, 36.0254} *)


      And of course, plotting shows that the eigenequation is not satisfied at all:



      With[{u = fl[[1]], b = λ[[1]]},
      Plot[(1 - ζ^2) u''[ζ] - ζ u'[ζ] + (b + 2 q - 4 q ζ^2) u[ζ], {ζ, 0, 1}]]
      With[{u = flt[[1]], b = λt[[1]]},
      Plot[(1 - ζ^2) u''[ζ] - ζ u'[ζ] + (b + 2 q - 4 q ζ^2) u[ζ], {ζ, 0, 1}]]


      What was wrong with my attempt? If I can get this example to work, I should be able to apply it to my actual, more complicated problem, so any Good Ideas would be welcome.










      share|improve this question









      New contributor




      宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
      Check out our Code of Conduct.







      $endgroup$




      To be able to apply the differentialequation capabilities of Mathematica to my graduate thesis, I am trying to apply NDEigensystem to an eigenproblem whose solution I know, but I am having some trouble doing so.



      As a test problem, I am using an algebraic version of the Mathieu equation,



      $$(1-zeta^{2})w^{primeprime}-zeta w^{prime}+left(a+2q-4qzeta^{2}right)w=0$$



      For this example I set $q=4/3$ and take only the first three eigenpairs:



      m = 3; q = 4/3;
      op = -(1 - ζ^2) u''[ζ] + ζ u'[ζ] + 2 q (2 ζ^2 - 1) u[ζ];
      bc = DirichletCondition[u[ζ] == 0, True];
      {λ, fl} = NDEigensystem[{op, bc}, u, {ζ, 0, 1}, m];


      I chose the Mathieu equation as a nontrivial example as Mathematica already has a function for it's evaluation:



      λt = Table[MathieuCharacteristicB[2 k, q], {k, m}];
      flt = Table[With[{j = j},
      MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];


      The problem is, I do not get the expected eigenvalues!



      λ
      (* {4.0708, 17.3259, 39.1877} *)
      N[λt]
      (* {3.85298, 16.0581, 36.0254} *)


      And of course, plotting shows that the eigenequation is not satisfied at all:



      With[{u = fl[[1]], b = λ[[1]]},
      Plot[(1 - ζ^2) u''[ζ] - ζ u'[ζ] + (b + 2 q - 4 q ζ^2) u[ζ], {ζ, 0, 1}]]
      With[{u = flt[[1]], b = λt[[1]]},
      Plot[(1 - ζ^2) u''[ζ] - ζ u'[ζ] + (b + 2 q - 4 q ζ^2) u[ζ], {ζ, 0, 1}]]


      What was wrong with my attempt? If I can get this example to work, I should be able to apply it to my actual, more complicated problem, so any Good Ideas would be welcome.







      differential-equations finite-element-method






      share|improve this question









      New contributor




      宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
      Check out our Code of Conduct.











      share|improve this question









      New contributor




      宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
      Check out our Code of Conduct.









      share|improve this question




      share|improve this question








      edited 1 hour ago









      user21

      21k55998




      21k55998






      New contributor




      宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
      Check out our Code of Conduct.









      asked 2 hours ago









      宮川園子宮川園子

      211




      211




      New contributor




      宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
      Check out our Code of Conduct.





      New contributor





      宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
      Check out our Code of Conduct.






      宮川園子 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
      Check out our Code of Conduct.






















          2 Answers
          2






          active

          oldest

          votes


















          3












          $begingroup$

          If you refine the mesh, you will get closer:



          m = 3; q = 4/3;
          op = -(1 - [Zeta]^2) u''[[Zeta]] + [Zeta] u'[[Zeta]] +
          2 q (2 [Zeta]^2 - 1) u[[Zeta]];
          bc = DirichletCondition[u[[Zeta]] == 0, True];
          {[Lambda], fl} =
          NDEigensystem[{op, bc}, u, {[Zeta], 0, 1}, m,
          Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions"
          -> {"MaxCellMeasure" -> 0.00001}}}}];

          [Lambda]
          {3.855, 16.074, 36.064}

          [Lambda]t = Table[MathieuCharacteristicB[2 k, q], {k, m}];
          flt = Table[
          With[{j = j},
          MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];

          [Lambda]t // N
          {3.852, 16.058, 36.025}





          share|improve this answer









          $endgroup$





















            2












            $begingroup$

            It looks to me like NDEigensystem is struggling with the singularity at $zeta=1$, as does the method that I'm going to show. But perhaps it'll be useful for you, at least as a cross-check.



            I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions, the example notebook on the github or this introduction for some more details.



            First we install the package (only need to do this the first time):



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


            Then we first need to turn the ODEs into a matrix form $mathbf{y}'=mathbf{A} cdot mathbf{y}$, using my function ToMatrixSystem:



            Needs["CompoundMatrixMethod`"]
            sys[ζend_] = ToMatrixSystem[op == a u[ζ], {u[0] == 0, u[ζend] == 0}, u, {ζ, 0, ζend}, a]


            Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $a$ and $zeta_{end}$; this is an analytic function whose roots coincide with eigenvalues of the original equation.



            Plugging in $zeta_{end} = 1$ fails due to the singularity, but you can try moving the endpoint slightly away:



            FindRoot[Evans[a, sys[1 - 10^-3]], {a, 3}]
            (* {a -> 4.00335} *)


            Moving the endpoint closer approaches the correct value, but I can't get the exact value with this method.



            FindRoot[Evans[a, sys[1 - 10^-10], WorkingPrecision -> 30], {a, 3}, 
            WorkingPrecision -> 30] // Quiet
            (* {a -> 3.85301} *)


            You can see the same effect for the other roots.






            share|improve this answer









            $endgroup$














              Your Answer








              StackExchange.ready(function() {
              var channelOptions = {
              tags: "".split(" "),
              id: "387"
              };
              initTagRenderer("".split(" "), "".split(" "), channelOptions);

              StackExchange.using("externalEditor", function() {
              // Have to fire editor after snippets, if snippets enabled
              if (StackExchange.settings.snippets.snippetsEnabled) {
              StackExchange.using("snippets", function() {
              createEditor();
              });
              }
              else {
              createEditor();
              }
              });

              function createEditor() {
              StackExchange.prepareEditor({
              heartbeatType: 'answer',
              autoActivateHeartbeat: false,
              convertImagesToLinks: false,
              noModals: true,
              showLowRepImageUploadWarning: true,
              reputationToPostImages: null,
              bindNavPrevention: true,
              postfix: "",
              imageUploader: {
              brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
              contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
              allowUrls: true
              },
              onDemand: true,
              discardSelector: ".discard-answer"
              ,immediatelyShowMarkdownHelp:true
              });


              }
              });






              宮川園子 is a new contributor. Be nice, and check out our Code of Conduct.










              draft saved

              draft discarded


















              StackExchange.ready(
              function () {
              StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f196891%2fusing-ndeigensystem-to-solve-the-mathieu-equation%23new-answer', 'question_page');
              }
              );

              Post as a guest















              Required, but never shown

























              2 Answers
              2






              active

              oldest

              votes








              2 Answers
              2






              active

              oldest

              votes









              active

              oldest

              votes






              active

              oldest

              votes









              3












              $begingroup$

              If you refine the mesh, you will get closer:



              m = 3; q = 4/3;
              op = -(1 - [Zeta]^2) u''[[Zeta]] + [Zeta] u'[[Zeta]] +
              2 q (2 [Zeta]^2 - 1) u[[Zeta]];
              bc = DirichletCondition[u[[Zeta]] == 0, True];
              {[Lambda], fl} =
              NDEigensystem[{op, bc}, u, {[Zeta], 0, 1}, m,
              Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions"
              -> {"MaxCellMeasure" -> 0.00001}}}}];

              [Lambda]
              {3.855, 16.074, 36.064}

              [Lambda]t = Table[MathieuCharacteristicB[2 k, q], {k, m}];
              flt = Table[
              With[{j = j},
              MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];

              [Lambda]t // N
              {3.852, 16.058, 36.025}





              share|improve this answer









              $endgroup$


















                3












                $begingroup$

                If you refine the mesh, you will get closer:



                m = 3; q = 4/3;
                op = -(1 - [Zeta]^2) u''[[Zeta]] + [Zeta] u'[[Zeta]] +
                2 q (2 [Zeta]^2 - 1) u[[Zeta]];
                bc = DirichletCondition[u[[Zeta]] == 0, True];
                {[Lambda], fl} =
                NDEigensystem[{op, bc}, u, {[Zeta], 0, 1}, m,
                Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions"
                -> {"MaxCellMeasure" -> 0.00001}}}}];

                [Lambda]
                {3.855, 16.074, 36.064}

                [Lambda]t = Table[MathieuCharacteristicB[2 k, q], {k, m}];
                flt = Table[
                With[{j = j},
                MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];

                [Lambda]t // N
                {3.852, 16.058, 36.025}





                share|improve this answer









                $endgroup$
















                  3












                  3








                  3





                  $begingroup$

                  If you refine the mesh, you will get closer:



                  m = 3; q = 4/3;
                  op = -(1 - [Zeta]^2) u''[[Zeta]] + [Zeta] u'[[Zeta]] +
                  2 q (2 [Zeta]^2 - 1) u[[Zeta]];
                  bc = DirichletCondition[u[[Zeta]] == 0, True];
                  {[Lambda], fl} =
                  NDEigensystem[{op, bc}, u, {[Zeta], 0, 1}, m,
                  Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions"
                  -> {"MaxCellMeasure" -> 0.00001}}}}];

                  [Lambda]
                  {3.855, 16.074, 36.064}

                  [Lambda]t = Table[MathieuCharacteristicB[2 k, q], {k, m}];
                  flt = Table[
                  With[{j = j},
                  MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];

                  [Lambda]t // N
                  {3.852, 16.058, 36.025}





                  share|improve this answer









                  $endgroup$



                  If you refine the mesh, you will get closer:



                  m = 3; q = 4/3;
                  op = -(1 - [Zeta]^2) u''[[Zeta]] + [Zeta] u'[[Zeta]] +
                  2 q (2 [Zeta]^2 - 1) u[[Zeta]];
                  bc = DirichletCondition[u[[Zeta]] == 0, True];
                  {[Lambda], fl} =
                  NDEigensystem[{op, bc}, u, {[Zeta], 0, 1}, m,
                  Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions"
                  -> {"MaxCellMeasure" -> 0.00001}}}}];

                  [Lambda]
                  {3.855, 16.074, 36.064}

                  [Lambda]t = Table[MathieuCharacteristicB[2 k, q], {k, m}];
                  flt = Table[
                  With[{j = j},
                  MathieuS[MathieuCharacteristicB[2 k, q], q, ArcCos[#]] &], {k, m}];

                  [Lambda]t // N
                  {3.852, 16.058, 36.025}






                  share|improve this answer












                  share|improve this answer



                  share|improve this answer










                  answered 57 mins ago









                  user21user21

                  21k55998




                  21k55998























                      2












                      $begingroup$

                      It looks to me like NDEigensystem is struggling with the singularity at $zeta=1$, as does the method that I'm going to show. But perhaps it'll be useful for you, at least as a cross-check.



                      I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions, the example notebook on the github or this introduction for some more details.



                      First we install the package (only need to do this the first time):



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


                      Then we first need to turn the ODEs into a matrix form $mathbf{y}'=mathbf{A} cdot mathbf{y}$, using my function ToMatrixSystem:



                      Needs["CompoundMatrixMethod`"]
                      sys[ζend_] = ToMatrixSystem[op == a u[ζ], {u[0] == 0, u[ζend] == 0}, u, {ζ, 0, ζend}, a]


                      Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $a$ and $zeta_{end}$; this is an analytic function whose roots coincide with eigenvalues of the original equation.



                      Plugging in $zeta_{end} = 1$ fails due to the singularity, but you can try moving the endpoint slightly away:



                      FindRoot[Evans[a, sys[1 - 10^-3]], {a, 3}]
                      (* {a -> 4.00335} *)


                      Moving the endpoint closer approaches the correct value, but I can't get the exact value with this method.



                      FindRoot[Evans[a, sys[1 - 10^-10], WorkingPrecision -> 30], {a, 3}, 
                      WorkingPrecision -> 30] // Quiet
                      (* {a -> 3.85301} *)


                      You can see the same effect for the other roots.






                      share|improve this answer









                      $endgroup$


















                        2












                        $begingroup$

                        It looks to me like NDEigensystem is struggling with the singularity at $zeta=1$, as does the method that I'm going to show. But perhaps it'll be useful for you, at least as a cross-check.



                        I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions, the example notebook on the github or this introduction for some more details.



                        First we install the package (only need to do this the first time):



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


                        Then we first need to turn the ODEs into a matrix form $mathbf{y}'=mathbf{A} cdot mathbf{y}$, using my function ToMatrixSystem:



                        Needs["CompoundMatrixMethod`"]
                        sys[ζend_] = ToMatrixSystem[op == a u[ζ], {u[0] == 0, u[ζend] == 0}, u, {ζ, 0, ζend}, a]


                        Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $a$ and $zeta_{end}$; this is an analytic function whose roots coincide with eigenvalues of the original equation.



                        Plugging in $zeta_{end} = 1$ fails due to the singularity, but you can try moving the endpoint slightly away:



                        FindRoot[Evans[a, sys[1 - 10^-3]], {a, 3}]
                        (* {a -> 4.00335} *)


                        Moving the endpoint closer approaches the correct value, but I can't get the exact value with this method.



                        FindRoot[Evans[a, sys[1 - 10^-10], WorkingPrecision -> 30], {a, 3}, 
                        WorkingPrecision -> 30] // Quiet
                        (* {a -> 3.85301} *)


                        You can see the same effect for the other roots.






                        share|improve this answer









                        $endgroup$
















                          2












                          2








                          2





                          $begingroup$

                          It looks to me like NDEigensystem is struggling with the singularity at $zeta=1$, as does the method that I'm going to show. But perhaps it'll be useful for you, at least as a cross-check.



                          I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions, the example notebook on the github or this introduction for some more details.



                          First we install the package (only need to do this the first time):



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


                          Then we first need to turn the ODEs into a matrix form $mathbf{y}'=mathbf{A} cdot mathbf{y}$, using my function ToMatrixSystem:



                          Needs["CompoundMatrixMethod`"]
                          sys[ζend_] = ToMatrixSystem[op == a u[ζ], {u[0] == 0, u[ζend] == 0}, u, {ζ, 0, ζend}, a]


                          Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $a$ and $zeta_{end}$; this is an analytic function whose roots coincide with eigenvalues of the original equation.



                          Plugging in $zeta_{end} = 1$ fails due to the singularity, but you can try moving the endpoint slightly away:



                          FindRoot[Evans[a, sys[1 - 10^-3]], {a, 3}]
                          (* {a -> 4.00335} *)


                          Moving the endpoint closer approaches the correct value, but I can't get the exact value with this method.



                          FindRoot[Evans[a, sys[1 - 10^-10], WorkingPrecision -> 30], {a, 3}, 
                          WorkingPrecision -> 30] // Quiet
                          (* {a -> 3.85301} *)


                          You can see the same effect for the other roots.






                          share|improve this answer









                          $endgroup$



                          It looks to me like NDEigensystem is struggling with the singularity at $zeta=1$, as does the method that I'm going to show. But perhaps it'll be useful for you, at least as a cross-check.



                          I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions, the example notebook on the github or this introduction for some more details.



                          First we install the package (only need to do this the first time):



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


                          Then we first need to turn the ODEs into a matrix form $mathbf{y}'=mathbf{A} cdot mathbf{y}$, using my function ToMatrixSystem:



                          Needs["CompoundMatrixMethod`"]
                          sys[ζend_] = ToMatrixSystem[op == a u[ζ], {u[0] == 0, u[ζend] == 0}, u, {ζ, 0, ζend}, a]


                          Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $a$ and $zeta_{end}$; this is an analytic function whose roots coincide with eigenvalues of the original equation.



                          Plugging in $zeta_{end} = 1$ fails due to the singularity, but you can try moving the endpoint slightly away:



                          FindRoot[Evans[a, sys[1 - 10^-3]], {a, 3}]
                          (* {a -> 4.00335} *)


                          Moving the endpoint closer approaches the correct value, but I can't get the exact value with this method.



                          FindRoot[Evans[a, sys[1 - 10^-10], WorkingPrecision -> 30], {a, 3}, 
                          WorkingPrecision -> 30] // Quiet
                          (* {a -> 3.85301} *)


                          You can see the same effect for the other roots.







                          share|improve this answer












                          share|improve this answer



                          share|improve this answer










                          answered 38 mins ago









                          KraZugKraZug

                          3,49821130




                          3,49821130






















                              宮川園子 is a new contributor. Be nice, and check out our Code of Conduct.










                              draft saved

                              draft discarded


















                              宮川園子 is a new contributor. Be nice, and check out our Code of Conduct.













                              宮川園子 is a new contributor. Be nice, and check out our Code of Conduct.












                              宮川園子 is a new contributor. Be nice, and check out our Code of Conduct.
















                              Thanks for contributing an answer to Mathematica Stack Exchange!


                              • Please be sure to answer the question. Provide details and share your research!

                              But avoid



                              • Asking for help, clarification, or responding to other answers.

                              • Making statements based on opinion; back them up with references or personal experience.


                              Use MathJax to format equations. MathJax reference.


                              To learn more, see our tips on writing great answers.




                              draft saved


                              draft discarded














                              StackExchange.ready(
                              function () {
                              StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f196891%2fusing-ndeigensystem-to-solve-the-mathieu-equation%23new-answer', 'question_page');
                              }
                              );

                              Post as a guest















                              Required, but never shown





















































                              Required, but never shown














                              Required, but never shown












                              Required, but never shown







                              Required, but never shown

































                              Required, but never shown














                              Required, but never shown












                              Required, but never shown







                              Required, but never shown







                              Popular posts from this blog

                              El tren de la libertad Índice Antecedentes "Porque yo decido" Desarrollo de la...

                              Castillo d'Acher Características Menú de navegación

                              Puerta de Hutt Referencias Enlaces externos Menú de navegación15°58′00″S 5°42′00″O /...