{VERSION 3 0 "IBM INTEL NT" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{PSTYLE " Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Text Output" -1 2 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 0 0 0 0 0 1 3 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 6 6 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading \+ 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 4 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Warning" 2 7 1 {CSTYLE "" -1 -1 "" 0 1 0 0 255 1 0 0 0 0 0 0 1 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 12 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 } {PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }} {SECT 0 {PARA 18 "" 0 "" {TEXT -1 33 "Eigenstructure of the Loop Schem e" }}{PARA 19 "" 0 "" {TEXT -1 25 "Denis Zorin, January 1998" }}{PARA 0 "" 0 "" {TEXT -1 309 "In this worksheet we compute the eigenvalues a nd eigenvectors of the subdivision matrix of the Loop scheme. This a nalysis repeats the original derivation of Loop from his MS thesis. \+ We also compute eigenvalues and eigenvectors of the modified Loop sche me, which always has largest eigenvalue equal to 1/2." }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 9 "Utilities" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(linalg):" }}{PARA 7 "" 1 "" {TEXT -1 32 "Warning , new definition for norm" }}{PARA 7 "" 1 "" {TEXT -1 33 "Warning, new definition for trace" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "re adlib(C):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "readlib(optimi ze):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "readlib(procbody): \+ readlib(procmake):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "read( `subdivmatrix-util.mpl`):" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 57 "Lo op scheme subdivision matrix, eigenvalues, eigenvectors" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "assume( c, real); additionally( c > = -1 ); additionally( c <= 1); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 301 "Loop := matrix( [[ (1 - K*alpha)*d, K*alpha* d, 0,0], \n [ (3/8)*d, 3/8 + (1/4)*c, \+ 0,0],\n [ (1/16)*d, 5/8 + (1/8)*c, 1/16, (1 /16)*(1+conjugate(omega))],\n [ (1/8)*d, (3/8) *(1 + omega), 0, 1/8]\n ]);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%LoopG-%'MATRIXG6#7&7&*&,&\"\"\"F,*&%\"KGF,%&alphaGF, !\"\"F,%\"dGF,*(F.F,F/F,F1F,\"\"!F37&,$F1#\"\"$\"\"),&F6F,%#c|irG#F,\" \"%F3F37&,$F1#F,\"#;,&#\"\"&F8F,F:#F,F8F?,&F?F,-%*conjugateG6#%&omegaG F?7&,$F1FD,&F6F,FIF6F3FD" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Loopvar := \{ omega = c + I*s\};" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#> %(LoopvarG<#/%&omegaG,&%#c|irG\"\"\"*&%\"IGF*%\"sGF*F*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "evalues := eigenvals(subs( \{ d = 0 , op(Loopvar)\}, evalm(Loop)));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%( evaluesG6&\"\"!,&#\"\"$\"\")\"\"\"%#c|irG#F+\"\"%#F+\"#;#F+F*" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "EV := max(evalues);" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%#EVG,&#\"\"$\"\")\"\"\"%#c|irG#F)\" \"%" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "LoopZero := map( eva lc, subs( \{ d = 1, c = 1, omega = 1\}, evalm(Loop))):" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 21 "eigenvals( LoopZero);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6&\"\"\"#F#\"\")#F#\"#;,&#\"\"&F%F#*&%\"KGF#%&alphaGF#!\" \"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 90 "Minimal value of the larges t eigenvalue of the first block is 1/4; determine he range for " } {XPPEDIT 18 0 "alpha" "6#%&alphaG" }{TEXT -1 2 "; " }{XPPEDIT 18 0 "al pha" "6#%&alphaG" }{TEXT -1 74 " clearly should be less than 5/8K an d greater than the following number:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "AlphaCrit := solve( 5/8 - K*alpha = 3/8 + 1/4*cos(2*Pi/K), alph a );\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%*AlphaCritG,$*&,&!\"\"\"\" \"-%$cosG6#,$*&%#PiGF)%\"KGF(\"\"#F)F)F0F(#F(\"\"%" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "evects := subs( s^2 = 1-c^2, [eigenvects(ma p(evalc,subs( \{ d = 0, op(Loopvar)\}, evalm(Loop))))]);" }}{PARA 12 " " 1 "" {XPPMATH 20 "6#>%'evectsG7&7%,&#\"\"$\"\")\"\"\"%#c|irG#F+\"\"% F+<#-%'VECTORG6#7&\"\"!F+,$*&,(F,\"#I\"#EF+*$F,\"\"#F.F+,(F:F.F,\"\"* \"\"&F+!\"\"#F+F;,$*(%\"IGF+,(*&FCF+F,F+F?%\"sGF+FCF?F+,&F,F+F+F+F?#F) F;7%#F+F*F+<#-F16#7&F4F4,$*&FCF+,(FCF+FEF+FFF+F+F?F+7%#F+\"#;F+<#-F16# 7&F4F4F+F47%F4F+<#-F16#7&F+F4F4F4" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 114 "for i from 1 to vectdim(eve cts) do \n if op(1,op(i, evects)) = EV then v := op(1, op(3, op(i, ev ects))); fi; \nod;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "Loop Eigenvector := map(simplify, map(evalc, eval(v)));" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#>%0LoopEigenvectorG-%'VECTORG6#7&\"\"!\"\"\"*&,&%#c|i rG\"\"#\"#8F*F*,&\"\"&F*F-\"\"%!\"\",$*&,(F-F*F*F**&%\"IGF*%\"sGF*F*F* ,&F-F*F*F*F3#\"\"$F." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 22 "Drop the \+ first element" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "LoopEigenvector := vector( [ seq(LoopEigenvector[i], i = 2..4)]); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%0LoopEigenvectorG-%'VECTORG6#7%\"\"\"*&,&%#c|irG\"\"# \"#8F)F),&\"\"&F)F,\"\"%!\"\",$*&,(F,F)F)F)*&%\"IGF)%\"sGF)F)F),&F,F)F )F)F2#\"\"$F-" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 13 "Generate code" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "OutputFile := `loop.cpp`: \n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "MakeClassHeader( Outp utFile, `Loop`, 1,2,2, TriBox):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 105 "Although eigenvalues and eigenvectors are very simple, for unifor mity with other schemes we generate code" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "fprintf( OutputFile,`virtual Float Eigenvalue( int K ) \{\\n return FR(3,8) + FR(1,4)*cos(Float(2)*Fpi()/Float(K)); \\n\} \\n\\n`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "GenerateEigenv ectorCode(LoopEigenvector, OutputFile):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 195 "fprintf( OutputFile,`virtual void EigenvalueRange( F loat c, Float& lambdamin, Float& lambdamax ) \{\\n`); \nfprintf( Outpu tFile,`\\n lambdamin = FR(3,8) + FR(1,4)*c; lambdamax = FR(5,8); \\n \}\\n\\n`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "fprintf( Out putFile,`\};`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "fclose(O utputFile):" }}}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 66 "Modified Loop s cheme subdivision matrix, eigenvalues, eigenvectors" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 101 "This is the exp ression for all blocks except 0th; the 0th block is the same as for th e standard Loop." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 267 "ModLoop := mat rix( [[ (1 - K*alpha)*d, K*alpha*d, 0,0],\n \+ [ (3/8)*d, (1/2)^m, 0,0],\n [ \+ (1/16)*d, 5/8 + (1/8)*c, 1/16, (1/16)*(1+conjugate(omega))],\n[ \+ (1/8)*d, (3/8)*(1 + omega), 0, 1/8]\n]);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(ModLoopG-%'MATRIXG6#7&7&*&,&\"\"\"F,*&%\"KGF,%&alpha GF,!\"\"F,%\"dGF,*(F.F,F/F,F1F,\"\"!F37&,$F1#\"\"$\"\"))#F,\"\"#%\"mGF 3F37&,$F1#F,\"#;,&#\"\"&F8F,%#c|irG#F,F8F?,&F?F,-%*conjugateG6#%&omega GF?7&,$F1FE,&F6F,FJF6F3FE" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "mevalues := eigenvals(subs( \{ d = 0, op(Loopvar)\}, evalm(ModLoop )));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%)mevaluesG6&\"\"!)#\"\"\"\" \"#%\"mG#F)\"#;#F)\"\")" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 24 "The ad missable range of " }{XPPEDIT 18 0 "alpha" "6#%&alphaG" }{TEXT -1 13 " is obvious. " }{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "mevects := subs( s^2 = 1-c^2, [eigenvects(map(evalc , subs( \{ d = 0, op(Loopvar) , m = 1\}, evalm(ModLoop))))]);" }} {PARA 11 "" 0 "" {TEXT -1 0 "" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%(me vectsG7&7%\"\"!\"\"\"<#-%'VECTORG6#7&F(F'F'F'7%#F(\"#;F(<#-F+6#7&F'F'F (F'7%#F(\"\"#F(<#-F+6#7&F'F(,&#\"#7\"\"(F(%#c|irG#\"\"%F?*&%\"IGF(,(*& FDF(F@F(!\"\"%\"sGF(FDFGF(7%#F(\"\")F(<#-F+6#7&F'F',$*&FDF(,(FDF(FFF(F HF(F(FGF(" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 119 "for i from 1 to vectdim(mevects) do \n if op(1,op (i, mevects)) = 1/2 then mv := op(1, op(3, op(i, mevects))); fi; \nod ;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "ModLoopEigenvector := map( eva lc, map(simplify, mv));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%3ModLoopE igenvectorG-%'VECTORG6#7&\"\"!\"\"\",&#\"#7\"\"(F*%#c|irG#\"\"%F.,(F/F *F*F**&%\"IGF*%\"sGF*F*" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 22 "Drop t he first element" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "ModLoopEigenvec tor := vector( [ seq(ModLoopEigenvector[i], i = 2..4)]); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%3ModLoopEigenvectorG-%'VECTORG6#7%\"\"\",&#\" #7\"\"(F)%#c|irG#\"\"%F-,(F.F)F)F)*&%\"IGF)%\"sGF)F)" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 13 "Generate code" }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 29 "OutputFile := `modloop.cpp`:\n" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 55 "MakeClassHeader( OutputFile, `ModLoop`, 2,4, 2, TriBox):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 105 "Although eigenval ues and eigenvectors are very simple, for uniformity with other scheme s we generate code" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 86 "fprintf( Outp utFile,`virtual Float Eigenvalue( int K ) \{\\n return FR(1,2); \\n\} \\n\\n`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 160 "fprintf( Outp utFile,`virtual void EigenvalueRange( Float c, Float& lambdamin, Float & lambdamax) \{\\n return lambdamin = FR(1,2); lambdamax = FR(1,2); \+ \\n\}\\n\\n`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "GenerateE igenvectorCode(ModLoopEigenvector, OutputFile):" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 26 "fprintf( OutputFile,`\};`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "fclose(OutputFile):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}}}}{MARK "5 8 0 0" 13 }{VIEWOPTS 1 1 0 3 2 1804 }