How to use NDEigenvalue to accurately estimate functional determinants?

Goal: Ultimately, I would like to find a trustworthy approximation for the ratio of the functional determinant of two differential operators using the formula $$\frac{\text{Det}[\hat{D}_0]}{\text{Det}[\hat{D}]} =\frac{\prod^\infty_i \lambda_{0,i}}{\prod^\infty_j \lambda_j} \approx \frac{\prod^N_i \lambda_{0,i}}{\prod^N_j \lambda_j}$$ where $$\lambda$$ are the eigenvalues of the respective operators and $$N$$ is the number of eigenvalues used in the approximation.

Step 1: defining the operators.

To define these two operators I first must solve the following coupled ODE boundary value problem in $$\phi(t)$$ and $$p_\phi(t)$$. The system depends on the parameter $$\tau > 0$$.

bvp = {D[p\[Phi][t],
t] == (-(1/\[Tau]))*
p\[Phi][t]^2*(Sin[2*\[Phi][t]]/2) + (1/\[Tau])*p\[Phi][t]*
Cos[2*\[Phi][t]] + (1/\[Tau])*(Sin[2*\[Phi][t]]/2),
D[\[Phi][t], t] ==
2*Pi - ((1/\[Tau])*Sin[2*\[Phi][t]])/2 + (1/\[Tau])*p\[Phi][t]*
Sin[\[Phi][t]]^2, \[Phi][0] ==
2*ArcTan[(-1 + Sqrt[1 + 4*Pi^2*\[Tau]^2])/(2*Pi*\[Tau])], \[Phi][
1] == 2*ArcTan[(-1 + Sqrt[1 + 4*Pi^2*\[Tau]^2])/(2*Pi*\[Tau])] +
2*Pi}


sol = ParametricNDSolveValue[bvp, {\[Phi][t], p\[Phi][t]}, {t, 0, 1}, {\[Tau]}]


Using these solutions I can define the two differential operators (which depend on a choice of parameter $$\tau$$ and the numerical solution of $$\phi(t)$$) with Dirichlet boundary conditions

$$\hat{D}_0$$

{-D[\[Tau]*Csc[sol[\[Tau]][[1]]]^2*D[\[Psi][t], t], t], \[Psi][0] ==
0, \[Psi][1] == 0}


and $$\hat{D}$$

{-D[\[Tau]*Csc[sol[\[Tau]][[1]]]^2*D[\[Psi][t], t], t] + (2*\[Tau]*Cot[sol[\[Tau]]
[[1]]]*Csc[sol[\[Tau]][[1]]]^2*D[sol[\[Tau]][[1]], {t, 2}] -
\[Tau]*Csc[sol[\[Tau]][[1]]]^4*(Cos[2*sol[\[Tau]][[1]]] + 2)*D[sol[\[Tau]][[1]], t]^2 +
2*Pi*Csc[sol[\[Tau]][[1]]]^4*(2*Pi*\[Tau]*(Cos[2*sol[\[Tau]][[1]]] + 2) - Sin[2*sol[\[Tau]]
[[1]]]))*\[Psi][t], \[Psi][0] == 0,
\[Psi][1] == 0}


Step 2: Numerically calculating the eigenvalues.

First selecting a particular value of the parameter $$\tau$$ to be $$0.1$$, I use the inbuilt function NDEigenvalue to numerically calculate N=39 Eigenvalues for $$\hat{D}_0$$ and multiply them together

Times @@ NDEigenvalues[{-D[0.1*Csc[sol[0.1][[1]]]^2*D[\[Psi][t], t], t], \[Psi][0] == 0,
\[Psi][1] == 0}, \[Psi][t], {t, 0, 1}, 39]


Similarly for $$\hat{D}$$

Times @@ NDEigenvalues[{-D[0.1*Csc[sol[0.1][[1]]]^2*D[\[Psi][t], t], t] +
(2*0.1*Cot[sol[0.1][[1]]]*Csc[sol[0.1][[1]]]^2*D[sol[0.1][[1]], {t, 2}] -
0.1*Csc[sol[0.1][[1]]]^4*(Cos[2*sol[0.1][[1]]] + 2)*D[sol[0.1][[1]], t]^2 +
2*Pi*Csc[sol[0.1][[1]]]^4*(2*Pi*0.1*(Cos[2*sol[0.1][[1]]] + 2) - Sin[2*sol[0.1][[1]]]))*\[Psi][t], \[Psi][0] == 0,
\[Psi][1] == 0}, \[Psi][t], {t, 0, 1}, 39]


Then I simply take the ratio. I do this also for additional choices of the parameter $$\tau$$.

The issue: The accuracy of this approximation is not good enough for small values of the parameter $$\tau$$.

From the underlying physics, I strongly expect that

• The ratio will be a monotonically increasing function of $$\tau$$

Unfortunately, after calculating a few of these ratios in the parameter range $$0<\tau<1$$, I find that the ratios do not appear to be monotonically increasing. Additionally, the order of magnitude seems to change significantly as I change $$N$$, the number of eigenvalues used in the approximation. These two things make me conclude the numerical solutions I have so far are untrustworthy.

If I try to increase the number of Eigenvalues used in the approximation beyond 39 I get an error message saying "A maximum number of 39 eigenvalues and functions can be computed for this discretized system."

Solution Attempt:

My main guess at how to solve this issue is to find a way to increase the number of eigenvalues that NDEigenvalue returns.

1. I have tried to increase the accuracy of the parametric solutions for $$\phi(t)$$, but this did not seem to help very much

2. I have also tried to change methods used in NDEigenvalues (see the code below)

NDEigenvalues[{-D[0.1*Csc[sol[0.1][[1]]]^2*D[\[Psi][t], t],t], \[Psi][0] == 0,
\[Psi][1] == 0}, \[Psi][t], {t, 0, 1}, 70, Method -> {"SpatialDiscretization" ->
{"FiniteElement","MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}]


This will return about 60 eigenvalues before returning another error message (see image below) but has the same issue as the original method and will also sometimes find eigenvalues that when compared with the set of the first 39 eigenvalues found by the default method are sometimes an order of magnitude different!

So far all of this code runs fairly quickly, so I feel that there should be a way to get Mathematica to think for longer and find more eigenvalues but so far I have not found how to do this in a satisfactory way. I do not know if there are other ways I should be trying to solve my problem and would very much appreciate any advice I can get on how to make progress.

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Jul 21 at 18:51
• There is a typo in the first line of you code with not open ( ) bvp = {D[p\[Phi][t], t] == (-(1/\[Tau]))*p\[Phi][t]^2*(Sin[2*\[Phi][t]]/2) + 1/\[Tau])*... Jul 21 at 22:56
• Thanks for letting me know, I will edit my question to fix this and other errors that may have crept in while formatting my question. Jul 21 at 23:05
• The typo issue should be resolved now! Jul 21 at 23:14
• What about using perturbation theory and a Jacobi formula i.e. det(A + hB) for small h is equal to det (A) + h Tr(adj(A)B) ? Jul 21 at 23:26