Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 72

Mathematica package for Lamé eigenvalues

BeginPackage["LameEigenvalues`"];

Unprotect[RealLameEigenvalueA, RealLameEigenvalueB];

ClearAll[RealLameEigenvalueA, RealLameEigenvalueB];

Begin["`Private`"];

RealLameEigenvalueA[(ν_Integer)?NonNegative, (n_Integer)?NonNegative, m_] /; n <= ν :=
Module[{res, M},
	res = Which[
		TrueQ[0 <= m <= 1],
			LameEigenvalueA[ν, n, m],
		TrueQ[m < 0],
			M = m/(m - 1);
			If[EvenQ[n],
				m*ν*(ν + 1) + (1 - m)*LameEigenvalueA[ν, n, M],
				m*ν*(ν + 1) + (1 - m)*LameEigenvalueB[ν, n, M]
			],
		TrueQ[m > 1],
			M = 1/m;
			If[EvenQ[ν - n],
				m*LameEigenvalueA[ν, n, M],
				m*LameEigenvalueB[ν, n + 1, M]
			],
		True,
			$Failed
	];
	res /; res =!= $Failed
];

RealLameEigenvalueB[(ν_Integer)?Positive, (n_Integer)?Positive, m_] /; n <= ν :=
Module[{res, M},
	res = Which[
		TrueQ[0 <= m <= 1],
			LameEigenvalueA[ν, n, m],
		TrueQ[m < 0],
			M = m/(m - 1);
			If[EvenQ[n],
				m*ν*(ν + 1) + (1 - m)*LameEigenvalueB[ν, n, M],
				m*ν*(ν + 1) + (1 - m)*LameEigenvalueA[ν, n, M]
			],
		TrueQ[m > 1],
			M = 1/m;
			If[EvenQ[ν - n],
				m*LameEigenvalueA[ν, n - 1, M],
				m*LameEigenvalueB[ν, n, M]
			],
		True,
			$Failed
	];
	res /; res =!= $Failed
];

End[];

SetAttributes[{RealLameEigenvalueA, RealLameEigenvalueB}, {Listable, NumericFunction, ReadProtected}];

Protect[RealLameEigenvalueA, RealLameEigenvalueB];

EndPackage[];

Last edited by lanxiyu (Today 00:22:21)