3
$\begingroup$

I define the following function which spits out a number S when I feed it a dxd square matrix:

energycoeffs[coeffmatrix_] := 
  Module[{psi, psinorm, U, P, dxU, dyU, L, S},
    psi[x_, y_] := coeffmatrix.thetabasis[x, y];
    psinorm[x_, y_] := Chop[Norm[psi[x, y]]];
    U[x_, y_] := 1/psinorm[x, y]* psi[x, y];
    P[x_, y_] := U[x, y].ConjugateTranspose[U[x, y]];
    dxU[x_, y_] := (U[x + dx/2, y] - U[x - dx/2, y])/dx /. dx -> 0.01 ;
    dyU[x_, y_] := (U[x, y + dy/2] - U[x, y - dy/2])/dy /. dy -> 0.01 ;
    L[x_, y_] := 
      Re[ConjugateTranspose[dxU[x, y]].dxU[x, y] + 
         ConjugateTranspose[dyU[x, y]].dyU[x, y] - 
         ConjugateTranspose[dxU[x, y]].P[x, y].dxU[x, y] - 
         ConjugateTranspose[dyU[x, y]].P[x, y].dyU[x, y]];
    S := 
      N1*N2*NIntegrate[
        L[x, y]*Boole[{x, y} ∈ ℛ], {x, 0, gam1 + Re[gam2]}, {y, 0, Im[gam2]}, 
        Method -> "AdaptiveMonteCarlo", AccuracyGoal -> 3];
  S]

I evaluate it a point and it works:

energycoeffs[{{1, 0, 0}, {0, 0, 0}, {0, 0, 1}}]
{{18.8666}}

But I try to use it to find the coefficients of the optimal diagonal matrix,

FindMinimum[
  {
    energycoeffs[{{p, 0, 0}, {0, q, 0}, {0, 0, r}}], 
    0 <= p <= 1.2, 0 <= q <= 1.2, 0 <= r <= 1.2
  }, 
  {p, q, r}]

it gives me the following error over and over again:

NIntegrate::inumr: The integrand Boole[{x, y} ∈ Polygon[{{0, 0}, {0, Power[<<2>>] π}, {Power[<<2>>] π, Power[<<2>>] π}, {Power[<<2>>] π, 0}}]] Re[<<1>>] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0, 1}, {0, 1}}. >>

Does anyone know what I am doing wrong here?

EDIT: Here is the full code leading up to the function.

(* Set parameters *)
d = 3; (* !! *)
N1 = 1; N2 = 1; (*Number of unit \
cells*)
tau = 
 Exp[I π/3]; (* γ2=tau*γ1*)
gam1 = π Sqrt[
  d]; gam2 = π Sqrt[d] tau;
thetap[x_, y_, p_] = (
  E^(-((I (Sqrt[d] π tau - 2 (x + I*y))^2)/(4 π tau)))
    EllipticTheta[
    3, -(((d + 2 p) π)/(2 d)) + (x + I*y)/(Sqrt[d] tau), 
    E^(-((I π)/(d tau)))])/Sqrt[-I d tau];
thetabasis[x_, y_] = Table[{thetap[x, y, p]}, {p, 0, d - 1}]; 
Dimensions[thetabasis[x, y]]


ℛ = 
  Polygon[{{0, 0}, {Re[gam2], Im[gam2]}, {gam1 + Re[gam2], 
     Im[gam2]}, {gam1, 0}}] ;
Graphics[ℛ]

UPDATE: c186282 hit the nail on the head. The output of the function was not a scalar but a list with the scalar as the [1,1] value. E.g.-

energycoeffs[{{1, 0, 0}, {0, 0, 0}, {0, 0, 1}}]
{{18.8666}}

So now, I just use Norm[] to extract the scalar value. However, now I have a further question. FindMinimum is no longer throwing up any error messages, but it is extremely slow! It's been running for 10 minutes now without giving an answer. Would anyone have any advice as to how to fix it? Is it a problem with how I define the function?

UPDATE: I tried out the solution suggested by @george2079, and it works pretty well. I created a discretised mesh.

nsamp = 14.; 
dx = (gam1 + Re[gam2])/nsamp  ; dy = 
 Im[gam2]/nsamp ;(* Equivalent dx value. Note: x and y directions are \
not same, depending on shape of the unit cell. *)

axis = Flatten[Table[{i*dx, j*dy}, {i, 0, nsamp}, {j, 0, nsamp}], 
  1] ;
Dimensions[axis]

(* Select points in unit cell, always <= L1 *)
pts = Select[axis, # \[Element] \[ScriptCapitalR] &];
Dimensions[pts]
ListPlot[pts]

Then I copy and pasted his newly defined Module[], with a slight tweak in the integration definition:

energycoeffs[coeffmatrix_] := 
 Module[{psi, psinorm, U, P, dxU, dyU, L, S},
  psi[x_, y_] := coeffmatrix.thetabasis[x, y];
  psinorm[x_, y_] := Chop[Norm[psi[x, y]]];
  U[x_, y_] := 1/psinorm[x, y]*psi[x, y];
  P[x_, y_] := U[x, y].ConjugateTranspose[U[x, y]];
  dxU[x_, y_] := 
   With[{dx = 0.01}, (U[x + dx/2, y] - U[x - dx/2, y])/dx];
  dyU[x_, y_] := 
   With[{dy = 0.01}, (U[x, y + dy/2] - U[x, y - dy/2])/dy];
  L[x_, y_] := 
   First@First@
     Re[ConjugateTranspose[dxU[x, y]].dxU[x, y] + 
       ConjugateTranspose[dyU[x, y]].dyU[x, y] - 
       ConjugateTranspose[dxU[x, y]].P[x, y].dxU[x, y] - 
       ConjugateTranspose[dyU[x, y]].P[x, y].dyU[x, y]];
  N1 N2 Total@(L @@@ pts) * dx *dy]

It works pretty well!

In[80]:= energycoeffs[IdentityMatrix[d]]  // Timing

Out[80]= {7.073925, 19.6902}

Interestingly, it takes 7x more time to get approximately the same answer that we got with randomly generated points!

$\endgroup$
11
  • 1
    $\begingroup$ can you include the definition for thetabasis, gam1, gam2, N1, and N2 so I could give it a try? $\endgroup$
    – Jason B.
    Feb 9, 2015 at 13:51
  • $\begingroup$ I think you should be looking at what happens when energycoeffs[{{p, 0, 0}, {0, q, 0}, {0, 0, r}}] is evaluated. $\endgroup$
    – m_goldberg
    Feb 9, 2015 at 14:04
  • $\begingroup$ If I scrape and paste your example directly then one problem I am seeing is Boole[{1, 2} \[Element] \[ScriptCapitalR]]. Mathematica does not understand that you dearly want Olde English reverse italic animated font with lots of cute curlies "R" to mean "Reals". When I replace that with Boole[{1, 2} \[Element] Reals] Mathematica understands you are asking if 1 and 2 are elements of the real numbers. I hope there is someone inside Wolfram who cries on a daily basis as a result of new typesetting features that make users think they MUST desktop publish EVERY detail of their Mathematica input. $\endgroup$
    – Bill
    Feb 9, 2015 at 19:47
  • $\begingroup$ Obviously I cannot run your code either but here are some ideas. Make sure that your function is returning a scalar and force your inputs to be numbers. energycoeffs[p_?NumericQ,q_?NumericQ,r_?NumericQ]:=energycoeffs[{{p, 0, 0}, {0, q, 0}, {0, 0, r}}] Then minimize energycoeffs[p,q,r] $\endgroup$
    – c186282
    Feb 10, 2015 at 1:02
  • $\begingroup$ I have edited the main post, but I have a side question as well. How exactly does the 'p_?NumericQ' syntax work? Is it some kind of typecasting? $\endgroup$
    – ap21
    Feb 11, 2015 at 19:47

1 Answer 1

2
$\begingroup$

Here is an approach that is much faster, but with potentially some loss of accuracy. Do away with NIntegrate and eval the function at preselected points in the region:

 nsamp = 200;
 pts = Select[ 
      Table[{RandomReal[{0, gam1 + Re[gam2]}], 
             RandomReal[{0, Im[gam2]}]}, {nsamp}] , # \[Element] \[ScriptCapitalR] & ];
 energycoeffs[coeffmatrix_] := 
    Module[{psi, psinorm, U, P, dxU, dyU, L, S},
        psi[x_, y_] := coeffmatrix.thetabasis[x, y];
        psinorm[x_, y_] := Chop[Norm[psi[x, y]]];
        U[x_, y_] := 1/psinorm[x, y]*psi[x, y];
        P[x_, y_] := U[x, y].ConjugateTranspose[U[x, y]];
        dxU[x_, y_] := With[{dx = 0.01}, (U[x + dx/2, y] - U[x - dx/2, y])/dx];
        dyU[x_, y_] := With[{dy = 0.01}, (U[x, y + dy/2] - U[x, y - dy/2])/dy];
        L[x_, y_] := First@First@
               Re[ConjugateTranspose[dxU[x, y]].dxU[x, y] + 
                  ConjugateTranspose[dyU[x, y]].dyU[x, y] -
                  ConjugateTranspose[dxU[x, y]].P[x, y].dxU[x, y] - 
                  ConjugateTranspose[dyU[x, y]].P[x, y].dyU[x, y]];
                   N1 N2 Total@(L @@@ pts) (gam1 + Re[gam2]) Im[gam2] / nsamp ]
        energycoeffs[{{1, 0, 0}, {0, 0, 0}, {0, 0, 1}}] // Timing

{1.263608, 19.7669}

compare with ,

{24.429757, 18.6672}

a 20 fold performance increase. I don't know how significant that error is for your application, but you might use this to get close to the minimum then use the integral form for refinement. ( The accuracy of this approach will be improved greatly if you can discretize (mesh) your region, instead of using random points )

$\endgroup$
0

Your Answer

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

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