\\=============== 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