(* ::Package:: *) BeginPackage["TrigInt`"] TrigInt::usage = "TrigInt is faster than Integrate for large trigonometric functions."; Equilateral::usage= "Eigenfunctions of Neumann and Dirichlet Laplacian on the equilateral triangle with vertices (0,0), (1,0), (1/2,Sqrt[3]/2): Equilateral[Dirichlet,Symmetric] [m,n] - 1<=m<=n Equilateral[Dirichlet,Antisymmetric] [m,n] - 1<=m(X psin[n][A,B,a,b]), X_. x^n_. Cos[A_.+B_. x]/;FreeQ[{X,A,B,n},x]:>(X pcos[n][A,B,a,b]), X_. Sin[A_.+B_. x]/;FreeQ[{X,A,B},x]->(X psin[0][A,B,a,b]), X_. Cos[A_.+B_. x]/;FreeQ[{X,A,B},x]->(X pcos[0][A,B,a,b]), X_. x^n_./;FreeQ[{X,n},x]:>(X pp[n][a,b]), X_/;FreeQ[{X},x]->X (b-a), X_:>(Message[TrigInt::nomatch,X];Integrate[X,{x,a,b}]) }; (* single integral *) TInt[f_,{x_,a_,b_}]:=( ff=f/.Sin[X_]:>Sin[Collect[X,{x},Simplify]]/.Cos[X_]:>Cos[Collect[X,{x},Simplify]]; ff=Expand[ff]; ff=ff+f0+f1; ff=Replace[ff,sub[x,a,b],1]; ff/.f0->0/.f1->0 ); (* eigenfunctions *) h=1; r=h/(2Sqrt[3]); u=r-Global`y; v=Sqrt[3]/2(Global`x-h/2)+(Global`y-r)/2; w=Sqrt[3]/2(h/2-Global`x)+(Global`y-r)/2; EqFun[f_,g_]:=f[Pi (-#1-#2)(u+2r)/(3r)]g[Pi (#1-#2)(v-w)/(9r)]+ f[Pi #1 (u+2r)/(3r)]g[Pi (2#2+#1)(v-w)/(9r)]+ f[Pi #2 (u+2r)/(3r)]g[Pi (-2#1-#2)(v-w)/(9r)]; Equilateral[Global`Neumann,Global`Symmetric]=Evaluate[Simplify[EqFun[Cos,Cos]]]&; Equilateral[Global`Neumann,Global`Antisymmetric]=Evaluate[Simplify[EqFun[Cos,Sin]]]&; Equilateral[Global`Dirichlet,Global`Symmetric]=Evaluate[Simplify[EqFun[Sin,Cos]]]&; Equilateral[Global`Dirichlet,Global`Antisymmetric]=Evaluate[Simplify[EqFun[Sin,Sin]]]&; Equilateral[Global`Eigenvalue]=Evaluate[4/27(Pi/r)^2(#1^2+#1 #2+#2^2)]&; Equilateral[Global`Neumann,Global`Half]=Equilateral[Global`Neumann,Global`Symmetric]/.Global`x->(1-Global`x)/2/.Global`y->Global`y/2//Simplify; Equilateral[Global`Dirichlet,Global`Half]=Equilateral[Global`Dirichlet,Global`Antisymmetric]/.Global`x->(1-Global`x)/2/.Global`y->Global`y/2//Simplify; Equilateral[Global`Eigenvalue,Global`Half]=Evaluate[1/27(Pi/r)^2(#1^2+#1 #2+#2^2)]&; Equilateral[]=T[1/2,Sqrt[3]/2]; Equilateral[Global`Half]=T[0,Sqrt[3]]; Square[Global`Dirichlet]=Evaluate[Sin[#1 Pi Global`x]Sin[#2 Pi Global`y]]&; Square[Global`Neumann]=Evaluate[Cos[#1 Pi Global`x]Cos[#2 Pi Global`y]]&; Square[Global`Eigenvalue]=Evaluate[Pi^2(#1^2+#2^2)]&; Square[Global`Dirichlet,Global`Half]=Evaluate[Sin[#1 Pi Global`x]Sin[#2 Pi Global`y]-(-1)^(#1+#2)Sin[#2 Pi Global`x]Sin[#1 Pi Global`y]]&; Square[Global`Neumann,Global`Half]=Evaluate[Cos[#1 Pi Global`x]Cos[#2 Pi Global`y]+(-1)^(#1+#2)Cos[#2 Pi Global`x]Cos[#1 Pi Global`y]]&; Square[Global`Half]=T[0,1]; (* Transplantation *) LT[{p1_,p2_,p3_},{q1_,q2_,q3_}]:=Module[{ff}, ff[x_]:={x.{aa,bb}+cc ,x.{dd,ee}+ff}; AffineTransform[{{{aa,bb},{dd,ee}},{cc,ff}}]/.Solve[{ff[p1]==q1,ff[p2]==q2,ff[p3]==q3},{aa,bb,cc,dd,ee,ff}][[1]] ]; ST[p_,q_]:=Thread[{Global`x,Global`y}->LT[p,q][{Global`x,Global`y}]]; Transplant[f_,T1_,T2_]:=f/.ST[T1,T2]; (* Rayleigh quotient *) Rayleigh[p__]:=Num[p]/Denom[p]; Num[f_,T__List]:=Total[GInt[f,#]&/@{T}]; Num[f_,Longest[T__List],p__]:=Num[f,T]+Num[p]; Denom[f_,T__List]:=Total[NInt[f,#]&/@{T}]; Denom[f_,Longest[T__List],p__]:=Denom[f,T]+Denom[p]; GInt[f_,T_]:=TrigInt[Del[f].Del[f],Limits[T]]; NInt[f_,T_]:=TrigInt[f^2,Limits[T]]; Limits[{{c_,0},{d_,0},{a_,b_}}]:=Sequence[{Global`y,Min[0,b],Max[0,b]},{Global`x,Min[c,d]+(a-Min[c,d])Global`y/b,Max[c,d]+(a-Max[c,d])Global`y/b}]; Limits[{{c_,0},{d_,0},{a_,b_},cond_}]:=Sequence@@Refine[{{Global`y,Min[0,b],Max[0,b]},{Global`x,Min[c,d]+(a-Min[c,d])Global`y/b,Max[c,d]+(a-Max[c,d])Global`y/b}},cond]; (* Limits[{{c_,0},{d_,0},{a_,b_}}]:=Sequence[{Global`y,0,b},{Global`x,c+(a-c)Global`y/b,d+(a-d)Global`y/b}]; *) Del[f_]:={D[f,Global`x],D[f,Global`y]}; Grad[f_]:=Del[f]; T[a_,b_]:={{0,0},{1,0},{a,b}}; T[a_,b_,c_]:={{0,0},{1,0},{a,b},c}; Area[{p1_,p2_,p3_}]:=Abs[Cross[Append[p2-p1,0],Append[p3-p1,0]]][[3]]/2//Simplify; Perimeter[l_]:=Total[Sqrt[#.#]&/@(l-RotateLeft[l])]//Simplify; End [ ] EndPackage [ ]