L-shape Eigenvalue - 100 Digits - Computer Program

Robert Stephen Jones
October, 2015
www.hbeLabs.com
rsjones7 at yahoo dot com
The following three files can be used to calculate the lowest Dirichlet eigenvalue of the Laplacian within the famous L-shape (two unit-edged squares joined to adjacent edges of a third) to 100 digits. The corresponding eigenfunction is used to create the familiar MathWorks logo (matlab).
The Helmholtz equation is, $$ \Delta\Psi({\mathbf r})+\lambda\Psi({\mathbf r}) = 0 $$ and the boundary condition is $\Psi({\mathbf r})=0$ on the L-shape boundary.
This program is designed to convince the reader that this technique works.

My computers are running Ubuntu linux, with pari/GP compiled with GMP and pthreads.
The first few lines of program detM.gp is where you need to customize for your computer.
In one terminal, run: "$ gp lsh.gp"
In another terminal, run "$ tailf results.dat" to watch the output.
There are comments in the routines.

lsh.gp

This is the main program:
\\=============== lsh.gp ================================ \\ $ gp lsh.gp # invoke and whatch the output file "results.dat" read("library.gp"); { \\ N=number of matching points. \\ Digits=estimate of precision based on known convergence rate comment = "<<< For N > 64, note the pure alternation, each providing a bound.>>>"; N = 4; \\ how low can we start? N must be even. rp = getrp(N); default(realprecision,rp); \\ this will be updated as required. Digits = N/2+15; \\ This will be refined. Guess of initial number of digits. tol = 10^(-(Digits-4)); \\ How precise to get the root first two roots. Ldn = 9.0; Lup = 10.0; \\ Initial bound to eigenvalue. We know it is near 9.64 Lvec=vector(0); \\ Store all intermediate roots as N increases. Lvec = concat(Lvec,mysolve(N,Lup,Ldn,tol)); \\ first root. N=N+2; Lvec = concat(Lvec,mysolve(N,Lup,Ldn,tol)); \\ second root. [Ldn,Lup]= findbounds(Lvec); \\ Need not actually be bounds. relerr = (Lup-Ldn)/(0.5*(Lup+Ldn)); \\ relative error (gap between these roots) Digits = round(-log(relerr)/log(10)); \\ Maybe, how many matching digits? tally=0; MatchingDigits=0; \\ until(relerr<10^-102, \\set your tolerance. e.g., "10^-52" for just over 50 digits until (MatchingDigits>=100, \\ ... or ... set number of matching digits N=N+2; \\ always increment N by 2, keep it even. demonstrates convergence rp = getrp(N); default(realprecision,rp); \\ set the realprecision based on N tol = 10^(-(max(1.2*Digits,Digits+10))); \\ find root to this much better than L Lup = precision(Lup,rp); Ldn = precision(Ldn,rp); Lvec = precision(Lvec,rp); \\ update precision of Lvec Lvec = concat(Lvec, mysolve(N,Ldn,Lup,tol)); \\ add the next root sgn1=if(Lvec[#Lvec]>Lvec[#Lvec-1],"+","-"); \\ does L increase or decrease sgn0=if(#Lvec>5,if(Lvec[#Lvec-1]>Lvec[#Lvec-2],"+","-"),""); \\ previous one if (sgn0=="+"&& sgn1=="-",msg="UPPER BOUND"; Lup=Lvec[#Lvec-1];tally++, sgn0=="-"&& sgn1=="+",msg="LOWER BOUND"; Ldn=Lvec[#Lvec-1];tally++, msg=""); if(tally<2, [Ldn,Lup]= findbounds(Lvec)); \\ find the best bounds so far relerr = (Lup-Ldn)/(0.5*(Lup+Ldn)); \\ update relative error Digits = (-log(relerr)/log(10)); \\ approximate correct digits if(tally>1,MatchingDigits=cntdgts(Ldn,Lup));\\ matching digits in the bound \\ most of the following is used to display the results in file:results.dat. digitsmsg=if(tally>1 && msg!="", Strprintf(" (%5.3g) %d",relerr,MatchingDigits),""); write("results.dat",Strprintf(" %s %s", msg, digitsmsg)); if(N==66,write("results.dat",comment)); fmt = Strprintf("%%3d %%%d.%df (%%s) ",Digits+9,Digits+7); write1("results.dat",Strprintf(fmt, N,Lvec[#Lvec],sgn1)); ); write("results.dat"); Lvup=Vec(Str(Lup));Lvdn=Vec(Str(Ldn)); comment = "The correct, unrounded %d digits of the lowest L-shape Dirichlet eigenvalue are:"; write("results.dat", Strprintf(comment,MatchingDigits)); write1("results.dat"," "); i=1;while(i<min(#Lvup,#Lvdn)&&Lvup[i]==Lvdn[i], write1("results.dat",Strprintf("%s",Lvup[i]));i++); write("results.dat"); quit(); }; \\=======================================================
[download]

detM.gp

This is a program that is spawned to calculate the point-matching matrix and its determinant:
\\=========== detM.gp =================================== \\ Customize for your computer here. This is dual core, 4GB RAM default(nbthreads,4); default(threadsize, 300 000 000); default(parisize, 3 000 000 000); \\------------------------------------------------------- read("library.gp"); { read("control.gp"); \\ read first time to use N rp = getrp(N); \\ make sure there is enough precision default(realprecision,rp); \\ setup precision read("control.gp"); \\ read again, to properly get L initialize(N); \\ set up m-values and matching point info M=getM(N,L); \\ get the point-matching matrix detM = matdet(M,1); \\ figure out its determinant fmtstr = Strprintf("%%%d.%df",rp+4,rp+1); write ("rawL.dat", Strprintf(fmtstr,L)); \\ store L fmtstr = Strprintf("%%%d.%dg",rp+3,rp); Dstr = Strprintf(fmtstr,detM); write ("rawD.dat", Dstr); \\ store the determinant quit(); \\ make sure this process terminates }; \\=======================================================
[download]

library.gp

This is a small, custom library file with a few subroutines:
\\============== very small library of routines used by lsh.gp and detM.gp ===== \\ Disclaimer: I'm certain there are better ways to make this happen.\ \\------------------------------------------------------------------------------ \\ Wrapper function for multithreading, note: globals in GP don't work here. WRAPFUN_ROWJ(N,j,m,krMP,sinMP)={parvector(N,i,besselj(m[i],krMP[j])*sinMP[i,j])} \\------------------------------------------------------------------------------ \\ Return the NxN matrix at this L value. rMP and sinMP are global getM(N,L)={ local(k=sqrt(L), krMP=k*rMP, M=matrix(N,N)); for(j=1,N,M[j,]=WRAPFUN_ROWJ(N,j,m,krMP,sinMP)); return(M);} \\------------------------------------------------------------------------------ \\ Spawn a GP subprocess to get the to calculate the matrix and determinant \\ file:control.gp is how N and L are passed to it. \\ file:rawD.dat is where results of the determinants are saved and values retrieved. getdet (N,L)={ local(fn="control.gp"); \\ temporary control file system(Strprintf("rm -rf %s",fn)); \\ possibly remove that temporary file write(fn,Strprintf("N = %d; L = %f;",N,L)); local(cmd=Strprintf("nice -n 19 gp -q detM.gp >> detM.log")); system(cmd); \\ run the linux command cmd local(D=readvec("rawD.dat")); \\ get the determinant value return(D[#D]);} \\ and return the most recent one added \\------------------------------------------------------------------------------ isodd(n)={!(n==2*floor(n/2))} \\ true if n is odd \\------------------------------------------------------------------------------ \\ This is a dumb but functional calculation of realprecision based on N getrp(N)={return(round(1.5*N))} \\------------------------------------------------------------------------------ \\ Set up the m-values and the matching points (MP) initialize(N)={ global(rMP,sinMP,m); \\ global except for wrapper function local(xMP,yMP,tMP); if (isodd(N),print("error, N=", N," must be even"); quit()); m = vector(N,i,(2/3)*(2*floor(3*i/2)-1)); \\ These are Chebyshev-like distributed matching points on two edges of a square \\ They are simple, yet provide very good (the best I have seen) convergence rates xMP = vector(N/2, i, 1); yMP = vector(N/2, i, (1/2)*(1-cos(Pi*i/(N/2)))); xMP = concat(xMP, vector(N/2, i, (1/2)*(1+cos(Pi*i/(N/2))))); yMP = concat(yMP, vector(N/2, i, 1)); rMP = vector(#xMP, i, sqrt(xMP[i]^2+yMP[i]^2)); tMP = vector(#xMP, i, if(xMP[i]>0, atan(yMP[i]/xMP[i]),xMP[i]==0,Pi/2)); sinMP = matrix(N,N,i,j,sin(m[i]*tMP[j])); kill(xMP);kill(yMP);kill(tMP); \\ does this really free up memory? return();} \\------------------------------------------------------------------------------ \\ Find the root of the NxN point-matching matrix determinant \\ bounded by L0 and L1 to a relative tolerance of tol. Each "getdet" \\ spawns a subprocess running "detM.gp". This routine uses simple \\ method to make sure root is in the interval, then a very simple \\ secant method, or, if not possible, splits interval in half. \\ Caution: Terseness makes this routine not very robust. mysolve(N,L0,L1,tol)={ L=[min(L0,L1),max(L0,L1),-1,-1]; D=vector(4); D[1] = getdet(N,L[1]); D[2] = getdet(N,L[2]); if(D[1]*D[2]>0, until(D[1]*D[2]<0, \\widen the gap if needed dL=L[2]-L[1]; L[1]=L[1]-dL; D[1]=getdet(N,L[1]); L[2]=L[2]+dL; D[2]=getdet(N,L[2]); ); ); until(dL < tol, \\ use a simple secant-based root-finding method slope = (D[2]-D[1])/(L[2]-L[1]); L[3] = L[1]-D[1]/slope; D[3] = getdet(N,L[3]); L[4] = L[3]-2*D[3]/slope; D[4] = getdet(N,L[4]); if(D[3]*D[4]<0, \\ use secant method if we can if(L[3]<L[4],L[1]=L[3];D[1]=D[3];L[2]=L[4];D[2]=D[4], L[1]=L[4];D[1]=D[4];L[2]=L[3];D[2]=D[3]) , \\ else split interval in half L[3]=(L[1]+L[2])/2; D[3] = getdet(N,L[3]); if(D[1]*D[3]<0, L[2]=L[3]; D[2]=D[3], L[1]=L[3]; D[1]=D[3]); ); \\ A challenge for those with even a little dyslexia, like me. dL=(L[2]-L[1])/L[1]; \\ relative error ); return ((L[1]+L[2])/2);} \\------------------------------------------------------------------------------ \\ return the most recent lower and upper bounds. Kludge that works only if \\ the convergence rate good, like for the L-shape. Works correctly if N > 68 \\ for the L-shape with Chebyshev-like matching points. The root-finding routine \\ (mysolve) will expand the interval a little so this mistake can be tolerated here. \\ However, for equally spaced matching points (per FHM-1967), it works from \\ very low (N=4) N-values, but the convergence rate is not as good. findbounds(Lvec)= { Lmin = min(Lvec[#Lvec],Lvec[#Lvec-1]); Lmax = max(Lvec[#Lvec],Lvec[#Lvec-1]); return ([Lmin,Lmax]); } \\------------------------------------------------------------------------------ \\ return the number of matching digits in the two numbers cntdgts(La, Lb)= { local(Va=Vec(Str(La)),Vb=Vec(Str(Lb)),i=1); while(i<=min(#Va,#Vb)&&(Va[i]==Vb[i]),i++); return(i-2); \\ one also for decimal point }
[download]

results.dat

This is what my computer output when run. The 50-digit value was reported within 2 minutes while it took another half an hour to reach 100 digits. See below for the meaning of the numbers.
8 9.63972074097 (-) 10 9.6397236366706 (+) 12 9.63972386144135 (+) UPPER BOUND 14 9.639723843845090 (-) LOWER BOUND (1.83 e-9) 8 16 9.639723843973879 (+) 18 9.639723844025192 (+) UPPER BOUND (1.87 e-11) 9 20 9.63972384402192146 (-) LOWER BOUND (3.39 e-13) 12 22 9.6397238440219306809 (+) 24 9.6397238440219417506 (+) UPPER BOUND (2.11 e-15) 14 26 9.639723844021941047890 (-) LOWER BOUND (7.29 e-17) 16 28 9.63972384402194105039216 (+) 30 9.63972384402194105287886 (+) UPPER BOUND (5.18 e-19) 17 32 9.6397238440219410527094777 (-) LOWER BOUND (1.76 e-20) 19 34 9.63972384402194105271094101 (+) 36 9.63972384402194105271150120 (+) UPPER BOUND (2.10 e-22) 20 38 9.6397238440219410527114585314 (-) LOWER BOUND (4.43 e-24) 22 40 9.639723844021941052711459142241 (+) 42 9.639723844021941052711459273606 (+) UPPER BOUND (7.70 e-26) 24 44 9.63972384402194105271145926207424 (-) LOWER BOUND (1.20 e-27) 26 46 9.639723844021941052711459262339596 (+) 48 9.639723844021941052711459262367698 (+) UPPER BOUND (3.04 e-29) 28 50 9.63972384402194105271145926236473386 (-) LOWER BOUND (3.08 e-31) 30 52 9.6397238440219410527114592623648168642 (+) 54 9.6397238440219410527114592623648240202 (+) UPPER BOUND (9.35 e-33) 31 56 9.639723844021941052711459262364823117685 (-) LOWER BOUND (9.36 e-35) 33 58 9.63972384402194105271145926236482315571561 (+) 60 9.63972384402194105271145926236482315643624 (+) UPPER BOUND (4.02 e-36) 35 62 9.639723844021941052711459262364823156260770 (-) LOWER BOUND (1.82 e-38) 37 64 9.63972384402194105271145926236482315626673875 (+) <<< For N > 64, note the pure alternation, each providing a bound.>>> 66 9.63972384402194105271145926236482315626738416 (+) UPPER BOUND (6.86 e-40) 39 68 9.6397238440219410527114592623648231562672825735 (-) LOWER BOUND (1.05 e-41) 40 70 9.63972384402194105271145926236482315626728985485 (+) UPPER BOUND (7.55 e-43) 42 72 9.6397238440219410527114592623648231562672895096840 (-) LOWER BOUND (3.58 e-44) 43 74 9.63972384402194105271145926236482315626728952771114 (+) UPPER BOUND (1.87 e-45) 44 76 9.639723844021941052711459262364823156267289525548882 (-) LOWER BOUND (2.24 e-46) 45 78 9.6397238440219410527114592623648231562672895258533445 (+) UPPER BOUND (3.16 e-47) 46 80 9.63972384402194105271145926236482315626728952581889523 (-) LOWER BOUND (3.57 e-48) 47 82 9.639723844021941052711459262364823156267289525822174357 (+) UPPER BOUND (3.40 e-49) 47 84 9.6397238440219410527114592623648231562672895258218819489 (-) LOWER BOUND (3.03 e-50) 48 86 9.63972384402194105271145926236482315626728952582190885081 (+) UPPER BOUND (2.79 e-51) 49 88 9.639723844021941052711459262364823156267289525821906213479 (-) LOWER BOUND (2.74 e-52) 51 90 9.6397238440219410527114592623648231562672895258219064806363 (+) UPPER BOUND (2.77 e-53) 52 92 9.63972384402194105271145926236482315626728952582190645366692 (-) LOWER BOUND (2.80 e-54) 53 94 9.639723844021941052711459262364823156267289525821906456350999 (+) UPPER BOUND (2.78 e-55) 54 96 9.6397238440219410527114592623648231562672895258219064560856402 (-) LOWER BOUND (2.75 e-56) 55 98 9.63972384402194105271145926236482315626728952582190645611197303 (+) UPPER BOUND (2.73 e-57) 55 100 9.639723844021941052711459262364823156267289525821906456109338840 (-) LOWER BOUND (2.73 e-58) 56 102 9.6397238440219410527114592623648231562672895258219064561096040182 (+) UPPER BOUND (2.75 e-59) 58 104 9.63972384402194105271145926236482315626728952582190645610957724207 (-) LOWER BOUND (2.78 e-60) 58 106 9.639723844021941052711459262364823156267289525821906456109579949476 (+) UPPER BOUND (2.81 e-61) 60 108 9.6397238440219410527114592623648231562672895258219064561095796753056 (-) LOWER BOUND (2.84 e-62) 61 110 9.63972384402194105271145926236482315626728952582190645610957970313419 (+) UPPER BOUND (2.89 e-63) 61 112 9.639723844021941052711459262364823156267289525821906456109579700301811 (-) LOWER BOUND (2.94 e-64) 63 114 9.6397238440219410527114592623648231562672895258219064561095797005908528 (+) UPPER BOUND (3.00 e-65) 64 116 9.63972384402194105271145926236482315626728952582190645610957970056128738 (-) LOWER BOUND (3.07 e-66) 65 118 9.639723844021941052711459262364823156267289525821906456109579700564317858 (+) UPPER BOUND (3.14 e-67) 66 120 9.6397238440219410527114592623648231562672895258219064561095797005640066112 (-) LOWER BOUND (3.23 e-68) 67 122 9.63972384402194105271145926236482315626728952582190645610957970056403864135 (+) UPPER BOUND (3.32 e-69) 68 124 9.639723844021941052711459262364823156267289525821906456109579700564035338662 (-) LOWER BOUND (3.43 e-70) 69 126 9.6397238440219410527114592623648231562672895258219064561095797005640356798601 (+) UPPER BOUND (3.54 e-71) 70 128 9.63972384402194105271145926236482315626728952582190645610957970056403564454640 (-) LOWER BOUND (3.66 e-72) 71 130 9.639723844021941052711459262364823156267289525821906456109579700564035648207815 (+) UPPER BOUND (3.80 e-73) 72 132 9.6397238440219410527114592623648231562672895258219064561095797005640356478275427 (-) LOWER BOUND (3.94 e-74) 72 134 9.63972384402194105271145926236482315626728952582190645610957970056403564786710306 (+) UPPER BOUND (4.10 e-75) 74 136 9.639723844021941052711459262364823156267289525821906456109579700564035647862980895 (-) LOWER BOUND (4.28 e-76) 75 138 9.6397238440219410527114592623648231562672895258219064561095797005640356478634110956 (+) UPPER BOUND (4.46 e-77) 75 140 9.63972384402194105271145926236482315626728952582190645610957970056403564786336613050 (-) LOWER BOUND (4.66 e-78) 76 142 9.639723844021941052711459262364823156267289525821906456109579700564035647863370837244 (+) UPPER BOUND (4.88 e-79) 77 144 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703438569 (-) LOWER BOUND (5.12 e-80) 79 146 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039564851 (+) UPPER BOUND (5.37 e-81) 80 148 9.639723844021941052711459262364823156267289525821906456109579700564035647863370390204503 (-) LOWER BOUND (5.65 e-82) 81 150 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703907774962 (+) UPPER BOUND (5.94 e-83) 82 152 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039071711008 (-) LOWER BOUND (6.26 e-84) 83 154 9.639723844021941052711459262364823156267289525821906456109579700564035647863370390723481948 (+) UPPER BOUND (6.61 e-85) 83 156 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703907228087800 (-) LOWER BOUND (6.98 e-86) 84 158 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039072287998226 (+) UPPER BOUND (7.39 e-87) 86 160 9.639723844021941052711459262364823156267289525821906456109579700564035647863370390722872442368 (-) LOWER BOUND (7.82 e-88) 87 162 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703907228732416943 (+) UPPER BOUND (8.29 e-89) 87 164 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039072287315686233 (-) LOWER BOUND (8.80 e-90) 88 166 9.639723844021941052711459262364823156267289525821906456109579700564035647863370390722873165875119 (+) UPPER BOUND (9.35 e-91) 89 168 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703907228731649165752 (-) LOWER BOUND (9.94 e-92) 90 170 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039072287316501862373 (+) UPPER BOUND (1.06 e-92) 90 172 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039072287316500774863 (-) LOWER BOUND (1.13 e-93) 92 174 9.639723844021941052711459262364823156267289525821906456109579700564035647863370390722873165008908693 (+) UPPER BOUND (1.20 e-94) 93 176 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703907228731650087848305 (-) LOWER BOUND (1.28 e-95) 94 178 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039072287316500879806789 (+) UPPER BOUND (1.37 e-96) 95 180 9.639723844021941052711459262364823156267289525821906456109579700564035647863370390722873165008796651903 (-) LOWER BOUND (1.47 e-97) 96 182 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703907228731650087968035026 (+) UPPER BOUND (1.57 e-98) 97 184 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039072287316500879678725782 (-) LOWER BOUND (1.69 e-99) 97 186 9.639723844021941052711459262364823156267289525821906456109579700564035647863370390722873165008796789000019 (+) UPPER BOUND (1.81 e-100) 99 188 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703907228731650087967888130195 (-) LOWER BOUND (1.94 e-101) 99 190 9.63972384402194105271145926236482315626728952582190645610957970056403564786337039072287316500879678883310740 (+) UPPER BOUND (2.08 e-102) 101 192 9.639723844021941052711459262364823156267289525821906456109579700564035647863370390722873165008796788830947797 (-) The correct, unrounded 101 digits of the lowest L-shape Dirichlet eigenvalue are: 9.6397238440219410527114592623648231562672895258219064561095797005640356478633703907228731650087967888
[download]
Notes
Example line:
 32 9.6397238440219410527094777 (-)   LOWER BOUND  (1.76 e-20) 19
At N=32, the determinant root "9.639..." decreased "(-)" compared to the previous one. This entry is effectively a "LOWER BOUND". The relative error between the (previous) upper bound and this lower bound is $\epsilon=1.76\times 10^{-20}$. There are 19 matching digits in these upper and lower bounds, "9.639723844021941052" The N=32 means there are 16 matching points (Chebyshev-like, see subroutine initialize) along each of two edges of one of the three squares.
Observe that for N>64, the values actually alternate, each providing one side of the next bound. When equally-spaced matching points are used instead, it will begin alternating from the beginning, but the convergence rate is a bit slower. I have confirmed this alternation up to about N=800, for both types of matching point distributions, which gives about 300 (equal-space) and 400 (Chebyshev) eigenvalue digits .

The calculation with N=192 is done simply to demonstrate that indeed the previous one is an upper bound.

If you want more digits, and faster, calculate a few digits, figure out the convergence rate, and jump to a higher N-value, perhaps picking up fifty digits on each iteration (not one). Repeat until you get as many digits as your computer will permit.

If you want this to work for you, make sure there is enough precision in the intermediate calculations to overcome the well-known "ill-conditioned" nature of this problem. The above program provides more than enough precision, and you can probably get by with a little less.
When I invoke the GP calculator, this is displayed:
                            GP/PARI CALCULATOR Version 2.7.3 (released)
                    amd64 running linux (x86-64/GMP-6.0.0 kernel) 64-bit version
                 compiled: Mar 20 2015, gcc version 4.8.2 (Ubuntu 4.8.2-19ubuntu1) 
                                     threading engine: pthread
                           (readline v6.3 enabled, extended help enabled)

                               Copyright (C) 2000-2015 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY 
WHATSOEVER.

Type ? for help, \q to quit.
Type ?12 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 500000
? 
www.hbeLabs.com
rsjones7 at yahoo dot com