1
$\begingroup$

I'm trying to implement the answer for "Integrate only one variable of a 2D interpolating function" (https://mathematica.stackexchange.com/a/161962/73672) but for interpolating functions returned from NDEigenSystem it isn't working.

First I solved for my required eigenfunctions and stored them as "funs":

ClearAll["Global`*"];
ClearAll[vals, funs, schröd];
A = 0.025;
Subscript[V, 0] = 1;
d = 2;
schröd = -A*d^2 π*D[ψ[n, φ], {φ, 2}] + 
   A/(4 π) (φ*φ*ψ[n, φ] + 
      2 I*D[ψ[n, φ], {n, 1}] - 
      D[ψ[n, φ], {n, 2}]) - 
   Subscript[V, 
     0] ((Cos[2 π*d*n]) + Cos[φ] - 20) ψ[
     n, φ];
Subscript[n, min] = -1/2; Subscript[n, max] = 1/2; 
Subscript[φ, min] = -π; 
Subscript[φ, max] = π;
Ω = 
  Rectangle[{Subscript[n, min], 
    Subscript[φ, min]}, {Subscript[n, max], 
    Subscript[φ, max]}];

{vals, funs} = 
  NDEigensystem[{schröd, 
    PeriodicBoundaryCondition[ψ[n, φ], 
     Subscript[φ, min] <= φ <= 
       Subscript[φ, max] && n == Subscript[n, max], 
     FindGeometricTransform[{{Subscript[n, min], 
         Subscript[φ, min]}, {Subscript[n, min], 
         Subscript[φ, max]}}, {{Subscript[n, max], 
         Subscript[φ, min]}, {Subscript[n, max], 
         Subscript[φ, max]}}][[2]]], 
    PeriodicBoundaryCondition[Exp[I 2 π n]*ψ[n, φ],
      Subscript[n, min] <= n <= Subscript[n, max] && φ == 
       Subscript[φ, max], 
     FindGeometricTransform[{{Subscript[n, min], 
         Subscript[φ, min]}, {Subscript[n, max], 
         Subscript[φ, min]}}, {{Subscript[n, min], 
         Subscript[φ, max]}, {Subscript[n, max], 
         Subscript[φ, 
          max]}}][[2]]]}, ψ, {n, φ} ∈ 
    Rectangle[{Subscript[n, min], 
      Subscript[φ, min]}, {Subscript[n, max], 
      Subscript[φ, max]}], 2];
Plot3D[{Evaluate[Abs[{funs[[1]][n, φ]  }^2]]}, {n, 
  Subscript[n, min], Subscript[n, max]}, {φ, 
  Subscript[φ, min], Subscript[φ, max]}, 
 PlotRange -> All, 
 PlotLabel -> 
  "Eig  fun    1 ", AxesLabel -> Automatic]

The answer to the other question was as follows:

da = 
  Flatten[
    Table[
      {t, tau, N @ Sin[2 (t + 3 tau)] Exp[-2 t - tau]}, 
      {t, 0, 2, 2/100}, {tau, 0, 5, 5/100}],
    1];
f = Interpolation @ da;
{{x1, x2}, {y1, y2}} = f["Domain"];
intx = Integrate[f[x, y], x] /. x -> x2;
    
nintx[y_?NumericQ] := Module[{x}, NIntegrate[f[x, y], {x, x1, x2}]];
    
Plot[nintx[y], {y, y1, y2}, PlotRange -> All]

However, trying to implement this myself gives an error after the third line, I believe because it reads n2 as 0.5:

f = funs[[1]];
{{n1, n2}, {φ1, φ2}} = f["Domain"];
intn = Integrate[f[n, φ], n] /. n -> n2;

General::ivar: 0.5 is not a valid variable.
Integrate::ilim: Invalid integration variable or limit(s) in 0.5.

$\endgroup$
4
  • $\begingroup$ funs doesn't evaluate with your code. Couldn't find definition of function \[Psi] $\endgroup$ Aug 11 '20 at 15:55
  • $\begingroup$ Where does funs come from? Have you checked that n doesn't have the value .5? $\endgroup$
    – m_goldberg
    Aug 11 '20 at 16:28
  • $\begingroup$ I've made some edits to the code and it is returning funs now. @UlrichNeumann are you referring to the fact that I am returning $\psi$ rather than $\psi(n,\varphi)$? $\endgroup$
    – Anthony
    Aug 11 '20 at 16:50
  • $\begingroup$ @m_goldberg funs is just what I'm storing the eigenfunctions returned from NDEigenSystem in. n is just a variable, but it looks like n2 gets set to 0.5 because that's the domain of the eigenfunction $\endgroup$
    – Anthony
    Aug 11 '20 at 16:53
0
$\begingroup$

Integrate won't integrate an InterpolatingFunction over an ElementMesh, which is what NDEigenSystem returns:

f@"ElementMesh"
(*  NDSolve`FEM`ElementMesh[{{-0.5, 0.5}, {-3.14159, 
     3.14159}}, {NDSolve`FEM`QuadElement["<" 408 ">"]}]  *)

We can integrate with Simpson's rule, since the mesh consists of quadratic QuadElements which have a node at the midpoints of the edges of the rectangles and in this case are evenly spaced:

nn = f["Grid"][[All, 1]] // DeleteDuplicates // Sort;
pp = f["Grid"][[All, 2]] // DeleteDuplicates // Sort;

(* check even spacing *)
{dn} = nn // Differences // DeleteDuplicates
(*  {0.0625}  *)
simp = With[{y = Outer[f, nn, pp]},
   (Total@y[[{1, -1}]] + 2 Total@y[[3 ;; -3 ;; 2]] + 4 Total@y[[2 ;; -2 ;; 2]])*
    dn/3
   ];
intn = Interpolation@Transpose@{pp, simp}

Compare with NIntegrate:

nint = Table[
   NIntegrate[f[n, phi], 
    Evaluate@Flatten@{n, First@f@"Domain"}], {phi, 
    f["Grid"][[All, 2]] // DeleteDuplicates // Sort}];
(simp - nint)/nint // Abs // Max
(*  2.00292*10^-15  *)
$\endgroup$
0
$\begingroup$

Try

intn[phi_?NumericQ] := NIntegrate[f[n, phi], {n, n1, n2}]

Table[{phi, intn[phi]}, {phi,Subdivide[\[CurlyPhi]1, \[CurlyPhi]2, 5]}] 
(*{{-3.14159, -0.00205648 - 0.00607021 I}, {-1.88496, -0.0253336 -0.0747785 I}, {-0.628319, -0.142423 -0.420398 I}, {0.628319, -0.142423 -0.420397 I}, {1.88496, -0.0253317 -0.0747728 I}, {3.14159, -0.00201903 - 0.00595968 I}}*)

Result is complex!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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