{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 Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "2 D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" 0 21 "" 0 1 0 0 0 1 0 0 0 0 2 0 0 0 0 }{CSTYLE "" -1 256 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 259 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "Courier" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "Cou rier" 0 1 0 0 0 0 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 3 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "H eading 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 8 4 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 8 2 0 0 0 0 0 0 -1 0 } {PSTYLE "" 2 6 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 2 0 0 }0 0 0 -1 -1 -1 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 "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 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 47 "Analysis of Subdivision based on C^4 box spline" }}{PARA 19 "" 0 "" {TEXT -1 21 "Denis Zorin, may 2000 " }}{SECT 0 {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 "Warn ing, new definition for trace" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 85 "with(share):\nwith(plots):\nwith(intpak);\nwith(codegen):\nreadl ib(optimize):\nreadlib(C):" }}{PARA 6 "" 1 "" {TEXT -1 70 "See ?share \+ and ?share,contents for information about the share library" }}{PARA 6 "" 1 "" {TEXT -1 22 "Share Library: intpak" }}{PARA 6 "" 1 "" {TEXT -1 48 "Authors: Connell, Amanda E. and Corless, Robert." }} {PARA 6 "" 1 "" {TEXT -1 41 "Description: Interval Arithmetic Package " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7#%%initG" }}{PARA 7 "" 1 "" {TEXT -1 35 "Warning, new definition for fortran" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 92 "This file contains some functions for code generat ion, and a fix for interval arithmetics " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "currentdir(\"D:\\\\users\\\\dzorin\\\\maple\");" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "read `subdivmatrix-util.mpl`;" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#QTc:\\Program~Files\\Maple~V~Release~5 ~-~Server\\BIN.WNT6\"" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 49 "Matrix initialization and barycenter computations" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 90 "Insert zero entries into the matrix: lazy subdivision, in itialization for primal schemes. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 108 "upsample := proc(M,i,j) if irem(i-1,2) = 0 and irem(j-1,2) = 0 th en M[(i-1)/2+1,(j-1)/2+1]; else 0; fi; end;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%)upsampleGR6%%\"MG%\"iG%\"jG6\"F*F*@%3/-%%iremG6$,&9% \"\"\"!\"\"F3\"\"#\"\"!/-F/6$,&9&F3F4F3F5F6&9$6$,&F2#F3F5F@F3,&F;F@F@F 3F6F*F*F*" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 254 "Primal to dual step : compute barycenters; shrinks the matrix size by 1 vertically. specia l treatment of top left entry, which depends on the center, which is p rovided as the second argument.\nMissing entries m[i,0] on the left r eplaced by omega*m[1,i-1]" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 206 "prim al2dual := proc(M,g,i,j) \n if j <> 1 \n then M[i,j-1]+M[i+1,j-1 ]+M[i,j]+M[i+1,j]; \n elif i <> 1 then omega*(M[1,i-1]+M[1,i])+M[i, j]+M[i+1,j]; \n else g+omega*M[1,1]+M[1,1]+M[2,1];\n fi;\nend;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,primal2dualGR6&%\"MG%\"gG%\"iG%\" jG6\"F+F+@'09'\"\"\",*&9$6$9&,&F.F/!\"\"F/F/&F26$,&F4F/F/F/F5F/&F26$F4 F.F/&F26$F9F.F/0F4F/,(*&%&omegaGF/,&&F26$F/,&F4F/F6F/F/&F26$F/F4F/F/F/ F:F/F " 0 " " {MPLTEXT 1 0 188 "dual2primal := proc(M,i,j) \n local m; \n if i <> 1 then \n M[i-1,j]+M[i-1,j+1]+M[i,j]+M[i,j+1];\n else \+ \n (M[j,1]+M[j+1,1])*conjugate(omega) + M[i,j]+M[i,j+1];\n fi; \+ \nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 64 "Two dualization steps t aking a primal sector matrix to the next " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 354 "makeSectorMatrixPrimal := proc(M, g, steps) local m, size, Mdual, Mprimal;\n global dual2primal, primal2dual; \n Mprimal := M;\n size := rowdim(M);\n for m from 0 to steps-1 do \n Mdual := matrix(size-m-1,size-m-1,(i,j)->primal2dual(Mprimal,g,i,j));\n \+ Mprimal := matrix(size-m-1,size-m-2, (i,j)->dual2primal(Mdual,i,j));\n od;\n evalm(Mprimal); \nend:" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 37 "Subdivision matrix from sector matrix" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 110 "This assumes N by N-1 matrix, which is numbered from the upper left corner 9see below for order illustration)." }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 188 "renumberQuadPrimal := proc(i,j,size) local ni , nj, l, n; \n nj := j; ni := i-1;\n l := max(ni,nj); \n if( ni \+ <= nj ) then n := i-1; else n := i+(i-j)-2; fi;\n l*(l-1) + n +1; \+ \n end;\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%3renumberQuadPrimalGR6 %%\"iG%\"jG%%sizeG6&%#niG%#njG%\"lG%\"nG6\"F/C'>8%9%>8$,&9$\"\"\"!\"\" F8>8&-%$maxG6$F5F2@%1F5F2>8'F6>FB,(F7\"\"#F3F9!\"#F8,(*&F;F8,&F;F8F9F8 F8F8FBF8F8F8F/F/F/" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 28 "Illustratio n of reordaring: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "matrix( 7,6, ( i,j)-> renumberQuadPrimal(i,j,5));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# -%'matrixG6#7)7(\"\"\"\"\"$\"\"(\"#8\"#@\"#J7(\"\"#\"\"%\"\")\"#9\"#A \"#K7(\"\"'\"\"&\"\"*\"#:\"#B\"#L7(\"#7\"#6\"#5\"#;\"#C\"#M7(\"#?\"#> \"#=\"#<\"#D\"#N7(\"#I\"#H\"#G\"#F\"#E\"#O7(\"#U\"#T\"#S\"#R\"#Q\"#P" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 175 "Make DFT subdivision matrix ou t of a sector matrix, with order of entries given by renumber. This as sumes that the initial sector matrix before subdivision had entries m[ i,j] " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 432 "MakeSubdivMatrix := proc( M, renumber) \n local S, i,j,r,q; \n S := matrix(coldim(M)*rowdim(M) , coldim(M)*rowdim(M));\n for i from 0 to rowdim(M)-1 do \n for j from 0 to coldim(M)-1 do\n for r from 0 to rowdim(M)-1 do \n \+ for q from 0 to coldim(M)-1 do \n S[renumber(i+1, j+ 1,rowdim(M)),renumber(r+1,q+1,rowdim(M))] := coeff(M[i+1,j+1],m[r+1,q+ 1]);\n od; \n od;\n od;\n od;\n evalm(S);\nend: \n \+ " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 "Sligtly different version fo r omega=0" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 508 "MakeSubdivMatrixZero \+ := proc(M, renumber) \n local S, i,j,r,q; \n S := matrix(coldim(M)*r owdim(M), coldim(M)*rowdim(M)+1);\n for i from 0 to rowdim(M)-1 do \+ \n for j from 0 to coldim(M)-1 do\n S[renumber(i+1, j+1,rowd im(M)),1] := coeff(M[i+1,j+1],m[0,0]);\n for r from 0 to rowdim( M)-1 do \n for q from 0 to coldim(M)-1 do \n S[renu mber(i+1, j+1,rowdim(M)),renumber(r+1,q+1,rowdim(M))+1] := coeff(M[i+1 ,j+1],m[r+1,q+1]);\n od; \n od;\n od;\n od;\n evalm( S);\nend: \n" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 31 "\\Subdivision \+ matrix computation" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 166 "Initial mat rix size is N by N-1; after upsampling, it's size increases to 2*N-1 \+ by 2*N-2 , after each pair of dualization steps, it decreases by 1 in \+ each dimension." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 7 "N := 4;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"NG\"\"%" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "M0 := \+ matrix( N,N-1,(i,j) -> m[i,j] );" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# >%#M0G-%'matrixG6#7&7%&%\"mG6$\"\"\"F-&F+6$F-\"\"#&F+6$F-\"\"$7%&F+6$F 0F-&F+6$F0F0&F+6$F0F37%&F+6$F3F-&F+6$F3F0&F+6$F3F37%&F+6$\"\"%F-&F+6$F EF0&F+6$FEF3" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 9 "Upsample." }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "M0 := matrix(2*N-1,2*N-2, (i,j)->u psample(M0,i,j-1));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#M0G-%'matrix G6#7)7(\"\"!&%\"mG6$\"\"\"F.F*&F,6$F.\"\"#F*&F,6$F.\"\"$7(F*F*F*F*F*F* 7(F*&F,6$F1F.F*&F,6$F1F1F*&F,6$F1F4F57(F*&F,6$F4F.F*&F,6$F4F1F*&F,6$F4 F4F57(F*&F,6$\"\"%F.F*&F,6$FGF1F*&F,6$FGF4" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 83 "Two dualization steps result in midpoint subdivision. g i s irrelevant in this case." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 104 "SMmi dpoint := map( simplify, subs(conjugate(omega) = 1/omega, makeSector MatrixPrimal(M0, m[0,0],1 )));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%+S MmidpointG-%'matrixG6#7(7'*&,(&%\"mG6$\"\"!F/\"\"\"*&%&omegaGF0&F-6$F0 F0F0\"\"#*&F,F0F2\"\"\"F0F7F2!\"\",$F3\"\"%,&F3F5&F-6$F0F5F5,$F " 0 " " {MPLTEXT 1 0 325 "baryDiag := proc(M,g,i,j) \n if i > 1 and j > 1 \+ \n then M[i-1,j]+M[i+1,j]+M[i,j+1]+M[i,j-1]; \n elif i > 1 and j = 1 then M[i-1,j]+M[i+1,j]+M[i,j+1]+M[1,i-1]*omega; \n elif i = 1 \+ and j > 1 then M[1,j-1]+M[1,j+1]+M[j+1,1]*conjugate(omega)+M[i+1,j]; \+ \n else g+M[1,j+1]+M[j+1,1]*conjugate(omega)+M[i+1,j]; \n fi; \n end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 128 "First filter step; the c entral vertex is unchanged by midoint subdivision; scale by 4 as we us e unscaled coeffs for barycenters." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 103 "SMbox1 := matrix(row dim(SMmidpoint)-1, coldim(SMmidpoint)-1, (i,j)->baryDiag(SMmidpoint,4* m[0,0],i,j)):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 235 "The effect of t he previous step on the center is to replase it with the barycenter o f edge neighbors, which is zero for omega <> 0, and (m[1,1]+m[0,0])/2 \+ otherwise; again, scale by 4, but twice; d is 1 for zero freq. and 0 \+ otherwise." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 102 "SMbox2 := matrix(rowdim(SMbox1)-1, coldim( SMbox1)-1, (i,j)->baryDiag(SMbox1,8*d*(m[1,1]+m[0,0]),i,j));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%'SMbox2G-%'matrixG6#7&7%,6*&%\"dG\" \"\",&&%\"mG6$F-F-F-&F06$\"\"!F4F-F-\"\")*&,(F2F-*&%&omegaGF-F/F-\"\"# *&F2F-F9\"\"\"F-F*&,.F6F-F8\"\"%FCFJF/F:F2F-F;F-F-FEFF/\"#;*&,*F2F-F/F-F8F-FCF-F-FEF*&,*F8F-FCF-*&F9 FFCF>FQF-FRF-F -FEF&F06$F-FTF:*&,&FQF:FRF:F-FEF*& ,,F8F:FCFJFQF:FRFJ*&F9F<,&F/F:F?F:F-F-F-FEFFAF-FC\"#=FRF:FUF>F8F>F2F-F;F-,8F/FaoF?FaoFOF-FC\"#9FUFhoFRF-&F06$F TF:F-FfnF-&F06$F:FTF-F2F-F8F-7%,8F6F-F8F5FCFLF/FJF2F-F;F-FQF:FRF>F\\oF -FUF:*&,,F6F-F/FJF?F:FAF-FCF:F-F9FF?F>FCFfoFUFfoFRF>FioF>F[pF:F8F:7%,8F2F-F/F-F8F5FCFaoFQF 5FRFao*&F9FF8FJFgpF:FioF>FQFJF\\oF-,:F/F-F?F-FCFhoFUFhoFRF hoFioFhoFgpF-&F06$FJF:F-F[pF-&F06$FTFTF-F8F-FQF-" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 168 "S := subs( d = 0, map( expand, map(simplify, subs( omega = cos(a) +I*sin(a),subs( conjugate(omega) = 1/omega, Mak eSubdivMatrix(SMbox2, renumberQuadPrimal, false))))));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"SG-%'matrixG6#7.7.,&\"#;\"\"\"-%$cosG6#%\"aG\" \"),(*&%\"IGF,-%$sinGF/F,!\"'F-\"\"'F8F,\"\"#\"\"!F:F:F:F:F:F:F:F:7.,, *(F4\"\"\"F-F,F5F>F9*$)F-F9F>F9\"#8F,F3FAF-\"#:,&\"#9F,F-F9,(F,F,F-F,F 3F,F,F:F,F:F:F:F:F:F:7.,&\"#CF,F-\"\"%,(F3!\")F-F1F1F,F1F9F:,&F3!\"#F- F9F:F:F:F:F:F:7.,(\"#=F,F-F8F3F8,(F3FMF-F9FPF,F8F8F:F9F:F:F:F:F:F:7.,( F1F,F-F1F3F1FH,(F9F,F-F9F3F9F1F9F1F:F:F:F:F:F:7.,(F8F,F-FPF3FP,(F3F9FP F,F-F9,&F-F8F3F8F9F:F8F:F:F:F:F:F:7.FPF2FPF8F:,&F3F7F-F8F9F:F:F:F:F:7. ,(FDF,F-F,F3F,,(F3!\"\"F-F,FDF,FDFDF,,(F3FhnF-F,F,F,F,F,F:F:F:F:7.,(F8 F,F-F9F3F9FPF8FPF8F8F:F9F:F:F:F:7.FEFDFEFDFDFDF:F,F,F:F,F,7.,(F9F,F-F8 F3F8FPFXF8F8FPF:F:F:F:F:F97.,(F,F,F-FDF3FDFfn,&F-FDF3FDFEF,FD,&F-F,F3F ,F:F:F:F:F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 142 "subs( omega^2*conjugate(omega) = o mega, subs( omega*conjugate(omega) = 1, map( expand, MakeSubdivMatrix( SMbox2, renumberQuadPrimal, false))));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7.7.,*%\"dG\"\" )\"#;\"\"\"-%*conjugateG6#%&omegaG\"\"%F0F1,&F-\"\"'F3F,\"\"#\"\"!F5F5 F5F5F5F5F5F57.,*\"#9F,F-F,F0F8*$)F0F4\"\"\"F,,(F-F,F8F,F0F,,&F,F,F0F,F ,F5F,F5F5F5F5F5F57.,(\"#CF,F-F4F0F4,&F-F*F*F,F*F4F5,$F-F4F5F5F5F5F5F57 .,&\"#=F,F0F3,&F-F4FEF,F3F3F5F4F5F5F5F5F5F57.,&F*F,F0F*F@,&F4F,F0F4F*F 4F*F5F5F5F5F5F57.,&F3F,F0FE,&FEF,F0F4,$F0F3F4F5F3F5F5F5F5F5F57.FEF2FEF 3F5,$F-F3F4F5F5F5F5F57.,&F8F,F0F,,&F-F,F8F,F8F8F,,&F-F,F,F,F,F,F5F5F5F 57.,&F3F,F0F4FEF3FEF3F3F5F4F5F5F5F57.F=F8F=F8F8F8F5F,F,F5F,F,7.,&F4F,F 0F3FEFMF3F3FEF5F5F5F5F5F47.,&F,F,F0F8FQ,$F0F8F=F,F8F0F5F5F5F5F," }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 156 "SZero := map( expand, map(simplify, subs( \{d = 1 , omega = 1\},subs( conjugate(omega) = 1/omega, MakeSubdivMatrixZero( SMbox2, renumberQuadPrimal, true)))));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%&SZeroG-%'matrixG6#7.7/\"#=\"#K\"#7\"\"#\"\"!F.F.F.F.F.F.F.F.7 /\"#9\"#I\"#;F-\"\"\"F.F3F.F.F.F.F.F.7/\"\")\"#GF2F5F-F.F-F.F.F.F.F.F. 7/\"\"'\"#C\"#?F8F8F.F-F.F.F.F.F.F.7/F-F2F9\"\"%F5F-F5F.F.F.F.F.F.7/F8 F9F:F8F-F.F8F.F.F.F.F.F.7/F-F*F,F*F8F.F8F-F.F.F.F.F.7/F3\"#:F@F0F0F3F- F3F3F.F.F.F.7/F.F5F*F8F*F8F8F.F-F.F.F.F.7/F.F-F0F-F0F0F0F.F3F3F.F3F37/ F.F5F*F8F8F8F*F.F.F.F.F.F-7/F3F@F@F0F-F3F0F3F.F.F.F.F3" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "SZero := stackmatrix( matrix( [[24,32,8,0,0,0,0,0,0,0,0,0,0]]),SZe ro); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%&SZeroG-%'matrixG6#7/7/\"#C \"#K\"\")\"\"!F-F-F-F-F-F-F-F-F-7/\"#=F+\"#7\"\"#F-F-F-F-F-F-F-F-F-7/ \"#9\"#I\"#;F1\"\"\"F-F6F-F-F-F-F-F-7/F,\"#GF5F,F1F-F1F-F-F-F-F-F-7/\" \"'F*\"#?F:F:F-F1F-F-F-F-F-F-7/F1F5F*\"\"%F,F1F,F-F-F-F-F-F-7/F:F*F;F: F1F-F:F-F-F-F-F-F-7/F1F/F0F/F:F-F:F1F-F-F-F-F-7/F6\"#:FAF3F3F6F1F6F6F- F-F-F-7/F-F,F/F:F/F:F:F-F1F-F-F-F-7/F-F1F3F1F3F3F3F-F6F6F-F6F67/F-F,F/ F:F:F:F/F-F-F-F-F-F17/F6FAFAF3F1F6F3F6F-F-F-F-F6" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "map( x -> x/64, [eigenvals(SZero)]);" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#7/#\"\"\"\"\"%F%#F%\"#k#F%\"#KF'F)#F% \"#;F+F+\"\"!F-F-F-" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "unas sign('S');" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 176 "lEV := subs( \{ cos(a) = c, sin(a) = s, _t[1]= 1\},simplify( linsolve( submatrix( evalm( (1/64)*transpose(S) - lambda * diag(seq(1,i=1..12))), 2..12,1. .12), [seq(0,i=1..11)])));" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%$lEVG- %'vectorG6#7.,$*&,6\"#;\"\"\"%'lambdaG!%?6%\"cG\"\"**$)F.\"\"#\"\"\"\" &[7#*$)F0F4F5F-*$)F.\"\"$F5!'k#R\"*$)F.\"\"%F5\"'W@E*&F:F5F0F-!%#>)*&F .F-F0F5!$#>*&F3F5F0F5\"%!G\"F5,.!\"\"F-F.\"#KF2!$?$FC!\"%FE\"#kF9\"%C5 !\"\"#F-\"\"),$*&,>*(%\"IGF-F3F5%\"sGF-!%k;FA!%Wh*&FUF5FVF5!\"#FE\"%k; FC!$G\"F0F;F7F-F9FXF2FenF.Ffn*(FUF5F:F5FVF5\"%Wh*(FUF5F0F5FVF5FH*(FUF5 F.F5FVF5\"$G\"F4F-F5FGFN#FHF4F-,$*&,6F.!$C#F2FFFjn\"#'*FC!#kF0F-FPF-F7 F-FinFHFE\"$o(FT!$o(F5FGFN#F-F,\"\"!,$*&,:F0!\"*FYFPFC\"$C#FE!%!G\"*&F .F5F8F5!#K**FUF5F.F5FVF5F0F5FIF7FHF.FaoF2FdoFTFFFinF-FjnF`oF5FGFN#FHF, FfoFfoFfoFfoFfoFfo" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 28 "Eigenvalu es and eigenvectors" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "CharP oly := collect( simplify( charpoly(S, lambda)), lambda);" }}{PARA 12 " " 1 "" {XPPMATH 20 "6#>%)CharPolyG,6*$)%'lambdaG\"#7\"\"\"\"\"\"*&,&-% $cosG6#%\"aG!#5!#cF+F+)F(\"#6F*F+*&,(\"$&)*F+*$)F.\"\"#F*!\")F.\"$'>F+ )F(\"#5F*F+*&,(F9\"$s#!%%*yF+F.!%i8F+)F(\"\"*F*F+*&,(F9!%g>F.\"%sQ\"&O E$F+F+)F(\"\")F*F+*&,*F.FIF9\"%S')*$)F.\"\"$F*\"$7&!&SW(F+F+)F(\"\"(F* F+*&,*F9!&OZ#FQ!%sI\"&OV*F+F.!&3M\"F+)F(\"\"'F*F+*&,*F.\"&/6$!&_>'F+FQ \"%cmF9\"&?.%F+)F(\"\"&F*F+*&,*\"&%Q;F+FQ!%WhF.!&Cm#F9!&oF$F+)F(\"\"%F *F+*&,(FQ\"%[?F.\"%#>)F9\"&S-\"F+)F(FSF*F+" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "CharP oly := factor( subs( cos(a) = c, CharPoly));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%)CharPolyG**),&%'lambdaG\"\"\"!\"#F)\"\"#\"\"\"),&F(F )!\"\"F)F+F,)F(\"\"$F,,@*$)F(\"\"&F,F)*&)F(\"\"%F,%\"cGF)!#5*$F7F,!#]* $F0F,\"$s'*&F0F,)F9F+F,!\")*&F0F,F9F,\"$O\"*&)F(F+F,F@F,\"$C#*$FEF,!%+ K*&FEF,F9F,!$;%*&F(F)F@F,!$7&*&F(F,F9F,FLF(\"%'4%F9\"%[?*$F@F,\"%gD*$) F9F1F,\"$7&F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 255 "CharPolyM ainFactor := collect( subs( lambda = 64*lambda, lambda^5-50*lambda^4-1 0*lambda^4*c+672*lambda^3+136*lambda^3*c-8*lambda^3*c^2+224*lambda^2*c ^2-3200*lambda^2-416*lambda^2*c+4096*lambda-512*lambda*c^2-512*lambda* c+512*c^3+2048*c+2560*c^2), lambda);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%3CharPolyMainFactorG,2*$)%' lambdaG\"\"&\"\"\"\"+C=ut5*&,&!*+3')Q)\"\"\"%\"cG!*g@xn\"F/)F(\"\"%F*F /*&,(F0\")%e^c$\"*o2;w\"F/*$)F0\"\"#F*!(_r4#F/)F(\"\"$F*F/*&,(F8\"'/v \"*!)+s58F/F0!(ORq\"F/)F(F:F*F/*&,(F0!&oF$\"'W@EF/F8FFF/F(F/F/F0\"%[?F 8\"%gD*$)F0F=F*\"$7&" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "Cha rPolyMainFactor := collect( CharPolyMainFactor/64^5, lambda);" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#>%3CharPolyMainFactorG,2*$)%'lambdaG\" \"&\"\"\"\"\"\"*&,&#!#D\"#KF+%\"cG#!\"&F0F+)F(\"\"%F*F+*&,(F1#\"#<\"$7 &#\"#@\"$G\"F+*$)F1\"\"#F*#!\"\"F:F+)F(\"\"$F*F+*&,(F>#\"\"(\"%#>)#F/ \"%[?F+F1#!#8FIF+)F(F@F*F+*&,(F1#FB\"&oF$#F+\"%'4%F+F>FQF+F(F+F+F1#F+ \"')GC&F>#F)\"(_r4#*$)F1FDF*#F+FX" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "latex(CharPolyMainFactor);" }}{PARA 6 "" 1 "" {TEXT -1 66 "\{\\lambda\}^\{5\}+\\left (-\{\\frac \{25\}\{32\}\}-\{\\frac \{ 5\}\{32\}\}\\,c\\right )\{" }}{PARA 6 "" 1 "" {TEXT -1 70 "\\lambda\}^ \{4\}+\\left (\{\\frac \{17\}\{512\}\}\\,c+\{\\frac \{21\}\{128\}\}-\{ \\frac \{1\}\{" }}{PARA 6 "" 1 "" {TEXT -1 70 "512\}\}\\,\{c\}^\{2\}\\ right )\{\\lambda\}^\{3\}+\\left (\{\\frac \{7\}\{8192\}\}\\,\{c\}^\{2 \}-" }}{PARA 6 "" 1 "" {TEXT -1 70 "\{\\frac \{25\}\{2048\}\}-\{\\frac \{13\}\{8192\}\}\\,c\\right )\{\\lambda\}^\{2\}+\\left (-" }}{PARA 6 "" 1 "" {TEXT -1 67 "\{\\frac \{1\}\{32768\}\}\\,c+\{\\frac \{1\}\{409 6\}\}-\{\\frac \{1\}\{32768\}\}\\,\{c\}^\{2\}" }}{PARA 6 "" 1 "" {TEXT -1 70 "\\right )\\lambda+\{\\frac \{1\}\{524288\}\}\\,c+\{\\frac \{5\}\{2097152\}\}\\,\{c\}^\{2\}+\{" }}{PARA 6 "" 1 "" {TEXT -1 28 " \\frac \{1\}\{2097152\}\}\\,\{c\}^\{3\}" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 148 "NumPts := 100: EigenvalsList := seq( sort(map( abs, \+ [solve( subs( c = evalf( 2*(n+1e-10)/NumPts),CharPolyMainFactor))])), \+ n = -NumPts/2..NumPts/2): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 117 "EigenvalsPlotLists := seq( [ seq( [-1 + 2*(i-1)/NumPts, abs(op(j, op(i,[EigenvalsList])))], i = 1..NumPts)] ,j=1..5):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "display(seq(plot(op(i,[EigenvalsPlotList s]),color=black), i = 1..5),color=black, axesfont=[TIMES,ITALIC,10], l abels=[``,``]);" }}{PARA 13 "" 1 "" {GLPLOT2D 752 477 477 {PLOTDATA 2 "6+-%'CURVESG6$7`q7$$!\"\"\"\"!$\"1++++++#>)!#G7$$!1+++++++)*!#;$\"1++ +RFle6!#>7$$!1+++++++'*F1$\"1+++TL?!H#F47$$!1*************R*F1$\"1+++t =j$R$F47$$!1+++++++#*F1$\"1+++8O(yY%F47$$!1+++++++!*F1$\"1+++c$>=^&F47 $$!1+++++++))F1$\"1+++I^JClF47$$!1+++++++')F1$\"1+++:J;/vF47$$!1++++++ +%)F1$\"1+++V?7]%)F47$$!1+++++++#)F1$\"1+++#>34O*F47$$!1+++++++!)F1$\" 1++++(>N-\"!#=7$$!1+++++++yF1$\"1+++[C;26F\\o7$$!1+++++++wF1$\"1+++M*y o=\"F\\o7$$!1+++++++uF1$\"1+++\"eDDE\"F\\o7$$!1+++++++sF1$\"1+++Lc&RL \"F\\o7$$!1+++++++qF1$\"1+++;&>5S\"F\\o7$$!1+++++++oF1$\"1+++e_cj9F\\o 7$$!1+++++++mF1$\"1+++H*Q9_\"F\\o7$$!1+++++++kF1$\"1+++e^[u:F\\o7$$!1+ ++++++iF1$\"1+++:xaA;F\\o7$$!1+++++++gF1$\"1+++W,Zl;F\\o7$$!1+++++++eF 1$\"1+++Ek4.s9F\\o7$$!1+++++++CF1$\"1+++ \"zK/S\"F\\o7$$!1+++++++AF1$\"1+++'34;K\"F\\o7$$!1+++++++?F1$\"1+++<)G dB\"F\\o7$$!1+++++++=F1$\"1+++cO\"G9\"F\\o7$$!1+++++++;F1$\"1+++Q*)*G/ \"F\\o7$$!1+++++++9F1$\"1+++#[L.O*F47$$!1+++++++7F1$\"1+++!*\\yA#)F47$ $!1+++++++5F1$\"1******yR3YKF47$$ \"1+++++++gF^z$\"1+++:(G!e\\F47$$\"1+++++++!)F^z$\"1+++^!\\qs'F47$$\"1 +++++++5F1$\"1+++7\"=pZJtF\\o7$$\"1+++++++ mF1$\"1+++^p#4h(F\\o7$$\"1+++++++oF1$\"1,++[%*>#*yF\\o7$$\"1+++++++qF1 $\"1+++s%=_<)F\\o7$$\"1+++++++sF1$\"1*****zI6*f%)F\\o7$$\"1+++++++uF1$ \"1+++\"R3iu)F\\o7$$\"1+++++++wF1$\"1+++tJ/M!*F\\o7$$\"1+++++++yF1$\"1 +++4?NB$*F\\o7$$Fb\\lF1$\"1+++`S29'*F\\o7$$\"1+++++++#)F1$\"1+++U6:1** F\\o7$$\"1+++++++%)F1$\"1+++oF&*>5F^z7$$\"1+++++++')F1$\"1+++b]T\\5F^z 7$$\"1+++++++))F1$\"1+++)*o**y5F^z7$$\"1+++++++!*F1$\"1+++[Mp36F^z7$$ \"1+++++++#*F1$\"1+++r+]Q6F^z7$$\"1*************R*F1$\"1+++SBTo6F^z7$$ \"1+++++++'*F1$\"1+++DgU)>\"F^z7$$\"1+++++++)*F1$\"1+++(3P&G7F^z-%'COL OURG6&%$RGBGF*F*F*-F$6$7`q7$F($\"1+++(*****\\iF^z7$F/$\"1+++/=tfbF^z7$ F6$\"1+++qzG#H&F^z7$F;$\"1+++)HIW4&F^z7$F@$\"1+++'y8E$\\F^z7$FE$\"1+++ *)p&Qz%F^z7$FJ$\"1+++LD^rYF^z7$FO$\"1+++U0khXF^z7$FT$\"1+++r8%F^z7$Fho$\"1+++U*3E1%F^z7$F]p$\"1+++C=P(*RF^z7$Fbp$\"1+++d.& e$RF^z7$Fgp$\"1+++X$Gx(QF^z7$F\\q$\"1+++=@uAQF^z7$Faq$\"1+++8)p1x$F^z7 $Ffq$\"1+++uLK@PF^z7$F[r$\"1+++VAauOF^z7$F`r$\"1+++A')=IOF^z7$Fer$\"1+ ++7O9)e$F^z7$Fjr$\"1+++eTI[NF^z7$F_s$\"1+++Q0e5NF^z7$Fds$\"1+++BV*[Z$F ^z7$Fis$\"1+++xmJp,KF^z7$F`w$\"1+++[Jv'=$F^z7$Few$\"1+++plMtJF^z7$Fjw$\"1 +++_AXhJF^z7$F_x$\"1+++$*)\\5:$F^z7$Fdx$\"1+++j)>@9$F^z7$Fix$\"1+++*3V Y8$F^z7$F^y$\"1+++15gGJF^z7$Fcy$\"1+++Ia(R7$F^z7$Fhy$\"1+++z&[27$F^z7$ F]z$\"1+++XI!*=JF^z7$Fbz$\"1+++0=U=JF^z7$Ffz$\"1+++V#)G>JF^z7$Fjz$\"1+ ++%>'[@JF^z7$F*$\"1++++++DJF^z7$Fb[l$\"1+++5Y\")HJF^z7$Fg[l$\"1+++3c\" f8$F^z7$F\\\\l$\"1+++R$*GVJF^z7$Fa\\l$\"1+++,I#>:$F^z7$Ff\\l$\"1+++VZ! =;$F^z7$F[]l$\"1+++pP#H<$F^z7$F`]l$\"1+++b/F&=$F^z7$Fe]l$\"1+++hk$))>$ F^z7$Fj]l$\"1+++!*[h8KF^z7$F_^l$\"1+++m.gHKF^z7$Fc^l$\"1+++.#*yYKF^z7$ Fh^l$\"1+++]&z^E$F^z7$F]_l$\"1+++Q:x%G$F^z7$Fb_l$\"1+++Nuc0LF^z7$Fg_l$ \"1+++M=dFLF^z7$F\\`l$\"1+++K=z]LF^z7$Fa`l$\"1+++)GP_P$F^z7$Ff`l$\"1++ +'*4#4S$F^z7$F[al$\"1+++***eyU$F^z7$F`al$\"1+++x32cMF^z7$Fdal$\"1+++'4 !e&[$F^z7$Fial$\"1+++%R9k^$F^z7$F^bl$\"1+++Qig[NF^z7$Fcbl$\"1+++fL>#e$ F^z7$Fhbl$\"1+++;%>sh$F^z7$F]cl$\"1+++sYt`OF^z7$Fbcl$\"1+++Ipz\"p$F^z7 $Fgcl$\"1+++uDZJPF^z7$F\\dl$\"1+++my$Gx$F^z7$Fadl$\"1+++y/)f\"QF^z7$Fe dl$\"1+++O9+hQF^z7$Fjdl$\"1+++*[TF^z7$Fcfl$\"1+++r! p'yTF^z7$Fhfl$\"1++++q\\TUF^z7$F]gl$\"1+++*H!)zI%F^z7$Fbgl$\"1+++IifyV F^z7$Ffgl$\"1+++hs%RX%F^z7$F[hl$\"1+++,3\"[`%F^z7$F`hl$\"1+++y!=Ai%F^z 7$Fehl$\"1+++r`ey]F^z7$Fiil$\"1+++%yXKC&F^z7$F^jl$\"1+++yCXhaF^zFbjl-F$6$ 7`qFijl7$F/$\"1+++8(eD+(F^z7$F6$\"1,++([!*=L(F^z7$F;$\"1+++2zD\"f(F^z7 $F@$\"1*****\\$4B9yF^z7$FE$\"1*****>/;Q,)F^z7$FJ$\"1,++*3'o'>)F^z7$FO$ \"1+++\"[1oO)F^z7$FT$\"1+++MDwE&)F^z7$FY$\"1+++t5Ny')F^z7$Fhn$\"1+++O; )G#))F^z7$F^o$\"1+++KMMh*)F^z7$Fco$\"1+++wi]%4*F^z7$Fho$\"1+++\")G)HA* F^z7$F]p$\"1+++65FZ$*F^z7$Fbp$\"1+++&H\"yn%*F^z7$Fgp$\"1+++jy&[e*F^z7$ F\\q$\"1+++HAz)p*F^z7$Faq$\"1+++)RM)4)*F^z7$Ffq$\"1+++$>+#=**F^z7$F[r$ \"1+++OyS-5F17$F`r$\"1+++RMw75F17$Fer$\"1+++0:!H-\"F17$Fjr$\"1+++e]$G. \"F17$F_s$\"1+++.edU5F17$Fds$\"1+++XU8_5F17$Fis$\"1+++;*>:1\"F17$F^t$ \"1+++>9uq5F17$Fct$\"1+++hm!)z5F17$Fht$\"1+++DGs)3\"F17$F]u$\"1+++^l\\ (4\"F17$Fbu$\"1++++R816F17$Fgu$\"1+++J0k96F17$F\\v$\"1+++X;-B6F17$Fav$ \"1+++D@GJ6F17$Ffv$\"1+++SkUR6F17$F[w$\"1+++^)eu9\"F17$F`w$\"1+++sKQb6 F17$Few$\"1+++NM?j6F17$Fjw$\"1+++eF#4<\"F17$F_x$\"1+++[Xay6F17$Fdx$\"1 +++c=2'=\"F17$Fix$\"1+++]w]$>\"F17$F^y$\"1+++!fa3?\"F17$Fcy$\"1+++v`63 7F17$Fhy$\"1+++2CH:7F17$F]z$\"1+++E\")QA7F17$Fbz$\"1+++7YSH7F17$Ffz$\" 1+++,TMO7F17$Fjz$\"1+++6%3KC\"F17$F*$\"1+++0(***\\7F17$Fb[l$\"1+++,;7R 7F17$Fg[l$\"1+++D!G$G7F17$F\\\\l$\"1+++(f8w@\"F17$Fa\\l$\"1+++1N(p?\"F 17$Ff\\l$\"1+++bES'>\"F17$F[]l$\"1+++@j*e=\"F17$F`]l$\"1+++w'\\a<\"F17 $Fe]l$\"1+++m!e];\"F17$Fj]l$\"1+++1nra6F17$F_^l$\"1+++**4UW6F17$Fc^l$ \"1+++Ai;M6F17$Fh^l$\"1+++:x%R7\"F17$F]_l$\"1+++'pgP6\"F17$Fb_l$\"1+++ 1/g.6F17$Fg_l$\"1+++e>Y$4\"F17$F\\`l$\"1+++3/M$3\"F17$Fa`l$\"1+++91Bt5 F17$Ff`l$\"1+++ct7j5F17$F[al$\"1+++q^-`5F17$F`al$\"1+++9%=H/\"F17$Fdal $\"1+++V6!G.\"F17$Fial$\"1+++armA5F17$F^bl$\"1+++q)4D,\"F17$Fcbl$\"1++ +HBK-5F17$Fhbl$\"1+++T0(4#**F^z7$F]cl$\"1+++K1E=)*F^z7$Fbcl$\"1+++Jp+: (*F^z7$Fgcl$\"1+++6b66'*F^z7$F\\dl$\"1,++*G$[1&*F^z7$Fadl$\"1+++>q*4S* F^z7$Fedl$\"1+++i3`%H*F^z7$Fjdl$\"1+++@V%p=*F^z7$F_el$\"1*****\\sy!y!* F^z7$Fdel$\"1+++^Svn*)F^z7$Fiel$\"1,++yLwb))F^z7$F^fl$\"1+++7p'=u)F^z7 $Fcfl$\"1+++pIyD')F^z7$Fhfl$\"1+++#Rxr])F^z7$F]gl$\"1+++(HYcQ)F^z7$Fbg l$\"1,++vZpg#)F^z7$Ffgl$\"1*****z-.<8)F^z7$F[hl$\"1+++^r(y*zF^z7$F`hl$ \"1,++<)p\"eyF^z7$Fehl$\"1+++D$\\6r(F^z7$Fjhl$\"1+++`^wavF^z7$F_il$\"1 +++s;*eQ(F^z7$Fdil$\"1*****pX=$*>(F^z7$Fiil$\"1+++*e_])pF^z7$F^jl$\"1* ****zH7#=nF^zFbjl-F$6$7`q7$F($\"1+++,+++DF17$F/$\"1+++:%[AG#F17$F6$\"1 +++EE,)>#F17$F;$\"1+++k/ZN@F17$F@$\"1+++vr0%3#F17$FE$\"1+++@apR?F17$FJ $\"1+++37I+?F17$FO$\"1+++m2kk>F17$FT$\"1+++*z7>$>F17$FY$\"1+++sDc,>F17 $Fhn$\"1+++9u=t=F17$F^o$\"1+++_V[Y=F17$Fco$\"1+++Z'=7#=F17$Fho$\"1+++S U?(z\"F17$F]p$\"1+++d1HuBas;F17$F[r$\"1++ +0MAa;F17$F`r$\"1+++\"))\\kj\"F17$Fer$\"1+++_!z\">;F17$Fjr$\"1+++^JP-; F17$F_s$\"1+++h%)*fe\"F17$Fds$\"1+++L[-q:F17$Fis$\"1+++!3DWb\"F17$F^t$ \"1+++(pu\"R:F17$Fct$\"1+++g8DC:F17$Fht$\"1+++&zM'4:F17$F]u$\"1+++)R1` \\\"F17$Fbu$\"1+++q\"\\7[\"F17$Fgu$\"1+++&QZuY\"F17$F\\v$\"1+++ul)QX\" F17$Fav$\"1+++mKbS9F17$Ffv$\"1+++X]VF9F17$F[w$\"1+++J-_99F17$F`w$\"1++ +3!)z,9F17$Few$\"1+++l\"e#*Q\"F17$Fjw$\"1+++<7*oP\"F17$F_x$\"1+++H\")o k8F17$Fdx$\"1+++d/k_8F17$Fix$\"1+++4,uS8F17$F^y$\"1+++t&z*G8F17$Fcy$\" 1+++9:N<8F17$Fhy$\"1+++7\"\\eI\"F17$F]z$\"1+++xcY%H\"F17$Fbz$\"1+++g]> $G\"F17$Ffz$\"1+++y5.s7F17$Fjz$\"1+++9\"o4E\"F17$F*$\"1+++&H++D\"F17$F b[l$\"1+++-/sc7F17$Fg[l$\"1+++t4Pj7F17$F\\\\l$\"1+++KR&*p7F17$Fa\\l$\" 1+++j1Zw7F17$Ff\\l$\"1+++&*G#HG\"F17$F[]l$\"1+++')>J*G\"F17$F`]l$\"1++ +C%RcH\"F17$Fe]l$\"1+++2l!>I\"F17$Fj]l$\"1+++aY638F17$F_^l$\"1+++b]E98 F17$Fc^l$\"1+++*)*e.K\"F17$Fh^l$\"1+++)f(RE8F17$F]_l$\"1+++x?QK8F17$Fb _l$\"1+++)[8$Q8F17$Fg_l$\"1+++;H>W8F17$F\\`l$\"1+++b8-]8F17$Fa`l$\"1++ +W)*zb8F17$Ff`l$\"1+++.$H:O\"F17$F[al$\"1+++\"p5sO\"F17$F`al$\"1+++%)[ %GP\"F17$Fdal$\"1+++)yK%y8F17$Fial$\"1+++E_(RQ\"F17$F^bl$\"1+++NIZ*Q\" F17$Fcbl$\"1+++()p#\\R\"F17$Fhbl$\"1+++')yL+9F17$F]cl$\"1+++ckq09F17$F bcl$\"1+++WM.69F17$Fgcl$\"1+++U&>jT\"F17$F\\dl$\"1+++jac@9F17$Fadl$\"1 +++g=xE9F17$Fedl$\"1+++(RR>V\"F17$Fjdl$\"1+++!poqV\"F17$F_el$\"1+++\"Q g@W\"F17$Fdel$\"1+++]]@Z9F17$Fiel$\"1+++)HLAX\"F17$F^fl$\"1+++$o:sX\"F 17$Fcfl$\"1+++!yi@Y\"F17$Fhfl$\"1+++>^2n9F17$F]gl$\"1+++QK&>Z\"F17$Fbg l$\"1+++Swzw9F17$Ffgl$\"1+++\\)3;[\"F17$F[hl$\"1+++QtQ'[\"F17$F`hl$\"1 +++,O8\"\\\"F17$Fehl$\"1+++'4[e\\\"F17$Fjhl$\"1+++'HJ0]\"F17$F_il$\"1+ ++RO=0:F17$Fdil$\"1+++qb!)4:F17$Fiil$\"1+++1vR9:F17$F^jl$\"1+++z)f*=:F 1Fbjl-F$6$7`qFd`o7$F/$\"1+++!)QhTFF17$F6$\"1+++D$z(\\GF17$F;$\"1+++%3< j$HF17$F@$\"1+++z/z6IF17$FE$\"1+++daF!3$F17$FJ$\"1+++:Y&Q9$F17$FO$\"1+ ++p.w.KF17$FT$\"1+++\\EzgKF17$FY$\"1+++a[]:LF17$Fhn$\"1+++:$)HoLF17$F^ o$\"1+++3ZZ>MF17$Fco$\"1+++^uEpMF17$Fho$\"1+++@8'y^$F17$F]p$\"1+++.bSl NF17$Fbp$\"1+++\\B-7OF17$Fgp$\"1+++pO\"yl$F17$F\\q$\"1+++T^'Gq$F17$Faq $\"1+++y'\\su$F17$Ffq$\"1+++Z)H5z$F17$F[r$\"1+++M)fU$QF17$F`r$\"1+++`o )p(QF17$Fer$\"1+++kBD>RF17$Fjr$\"1+++UH4hRF17$F_s$\"1+++s5a-SF17$Fds$ \"1+++-diVSF17$Fis$\"1+++\"zsV3%F17$F^t$\"1+++gc![7%F17$Fct$\"1+++&RX \\;%F17$Fht$\"1+++y5\"[?%F17$F]u$\"1+++v+UWUF17$Fbu$\"1+++)>)y$G%F17$F gu$\"1+++G*HHK%F17$F\\v$\"1+++=&e=O%F17$Fav$\"1+++phe+WF17$Ffv$\"1+++' 4C\"RWF17$F[w$\"1+++wE[xWF17$F`w$\"1+++%[rc^%F17$Few$\"1+++:%*p`XF17$F jw$\"1+++4Zd\"f%F17$F_x$\"1+++X]IHYF17$Fdx$\"1+++kv*om%F17$Fix$\"1+++V *eVq%F17$F^y$\"1+++.apTZF17$Fcy$\"1+++$y7*yZF17$Fhy$\"1+++Sl,;[F17$F]z $\"1+++/=,`[F17$Fbz$\"1+++$Q.**)[F17$Ffz$\"1+++5epE\\F17$Fjz$\"1+++TLR j\\F17$F*Fhbl7$Fb[l$\"1+++o&>l.&F17$Fg[l$\"1+++Kc&H2&F17$F\\\\l$\"1+++ m:J4^F17$Fa\\l$\"1+++!e!fX^F17$Ff\\l$\"1+++%p&z\"=&F17$F[]l$\"1+++)yHz @&F17$F`]l$\"1+++!e&*RD&F17$Fe]l$\"1+++ic***G&F17$Fj]l$\"1+++%[KfK&F17 $F_^l$\"1+++)Q3=O&F17$Fc^l$\"1+++*eDwR&F17$Fh^l$\"1*****z?'QLaF17$F]_l $\"1+++_A4paF17$Fb_l$\"1+++\\cu/bF17$Fg_l$\"1+++D#[.a&F17$F\\`l$\"1+++ L?_'F17$Fehl$\"1+++C\"ymb'F17$Fjhl$\"1+++ " 0 "" {MPLTEXT 1 0 164 "EV := subs( \{ cos(a) = c, sin(a) = s, _t[1]= 1\},simplify( \+ linsolve( submatrix( evalm( (1/64)*S - lambda * diag(seq(1,i=1..12))), 2..12,1..12), [seq(0,i=1..11)]))):" }}{PARA 12 "" 1 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "EVRe := map (simplify, ma p (evalc, map (Re, EV ))):" }}{PARA 12 "" 1 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "EVIm := map (simplify, map( \+ x -> x/s, map (evalc, map (Im, EV )))):" }}{PARA 12 "" 1 "" {TEXT -1 0 "" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 43 "Eigenvalue analysis and \+ interval estimation" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 155 "Eigenvalues cannot be found explicitly (degree 5 equat ion). Here obtain information about the eigenvalues to verify C1-con tinuity. We prove that for any " }{XPPEDIT 18 0 "m, k" "6$%\"mG%\"kG" }{TEXT -1 3 ", " }{XPPEDIT 18 0 "m = 1..k-1" "6#/%\"mG;\"\"\",&%\"kG \"\"\"\"\"\"!\"\"" }{TEXT -1 60 " the largest eigenvalue is real and \+ unique, and that for " }{XPPEDIT 18 0 " m <> k-1,1" "6$0%\"mG,&%\"kG \"\"\"\"\"\"!\"\"\"\"\"" }{TEXT -1 76 " the largest eigenvalue is less than the largest eigenvalue of blocks 1 and " }{XPPEDIT 18 0 "k-1" "6 #,&%\"kG\"\"\"\"\"\"!\"\"" }{TEXT -1 119 ". We also show that the uni que largest eigenvalue is a single eigenvalue in the interval [0.47+0. 2c, 0.52+0.2c], for " }{XPPEDIT 18 0 "k" "6#%\"kG" }{TEXT -1 12 " > 3 , where " }{XPPEDIT 18 0 "c = cos(2*Pi/k);" "6#/%\"cG-%$cosG6#*(\"\"# \"\"\"%#PiGF*%\"kG!\"\"" }{TEXT -1 7 " \nFor " }{XPPEDIT 18 0 "k = 3 " "6#/%\"kG\"\"$" }{TEXT -1 157 ", eigenvalues are examined separately . The proof is performed in several steps: \n(1). We show that for \+ c < 0, all roots of the characteristic polynomial " }{XPPEDIT 18 0 "P (c,lambda)" "6#-%\"PG6$%\"cG%'lambdaG" }{TEXT -1 169 " are less than 0 .51 (actually, they are less than 0.5, but due to numerical nature of \+ our calculations, we have to relax the upper boundary). \n(2). We show that for any " }{XPPEDIT 18 0 "c = 0..1" "6#/%\"cG;\"\"!\"\"\"" } {TEXT -1 30 ", there is a unique real root " }{XPPEDIT 18 0 "mu" "6#%# muG" }{TEXT -1 57 " in the interval [0.47+0.2c, 0.5+0.2c], and the fun ction " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" }{TEXT -1 117 " is C1 -continuous and increases. \n(3). We \"deflate\" the characteristic \+ polynomial (that is, divide by the monomial " }{XPPEDIT 18 0 "lambda - mu" "6#,&%'lambdaG\"\"\"%#muG!\"\"" }{TEXT -1 26 ") in symbolic form , with " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 50 " as indeterminates. Next, we verify t hat for all " }{XPPEDIT 18 0 "c = 0 .. 1;" "6#/%\"cG;\"\"!\"\"\"" } {TEXT -1 6 ", and " }{XPPEDIT 18 0 "mu = .47+.2*c .. .52+.2*c;" "6#/%# muG;,&$\"#Z!\"#\"\"\"*&$\"\"#!\"\"F*%\"cGF*F*,&$\"#_!\"#F**&$\"\"#!\" \"F*F/F*F*" }{TEXT -1 148 " , all roots of the deflated polynomial ar e inside the circle of radius 0.5 centered at 0 in the complex plane, \+ that is, have magnitudes less than " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6 #%\"cG" }{TEXT -1 10 " for any " }{XPPEDIT 18 0 "c > 0" "6#2\"\"!%\"c G" }{TEXT -1 10 ". \nAs for " }{XPPEDIT 18 0 " k" "6#%\"kG" }{TEXT -1 6 " > 4, " }{XPPEDIT 18 0 "cos(2*Pi/k) > 0.51" "6#2$\"#^!\"#-%$cosG6#* (\"\"#\"\"\"%#PiGF,%\"kG!\"\"" }{TEXT -1 63 ", the largest eigenvalue \+ cannot possibly correspond to a block " }{XPPEDIT 18 0 "m" "6#%\"mG" } {TEXT -1 12 ", for which " }{XPPEDIT 18 0 "cos(2*m*Pi/k) <= 0\n" "6#1- %$cosG6#**\"\"#\"\"\"%\"mGF)%#PiGF)%\"kG!\"\"\"\"!" }{TEXT -1 72 " . \+ From (3), it follows that the largest root has to be the real root " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" }{TEXT -1 10 " for some " } {XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 14 ". As for any " }{XPPEDIT 18 0 "m > 1, m < k-1" "6$2\"\"\"%\"mG2F%,&%\"kG\"\"\"\"\"\"!\"\"" } {TEXT -1 2 ", " }{XPPEDIT 18 0 "cos(2*m*Pi/k) < cos(2*Pi/k)" "6#2-%$co sG6#**\"\"#\"\"\"%\"mGF)%#PiGF)%\"kG!\"\"-F%6#*(\"\"#F)F+F)F,F-" } {TEXT -1 30 " , and we have shown (1) that " }{XPPEDIT 18 0 "mu(c)" "6 #-%#muG6#%\"cG" }{TEXT -1 25 ", increases, and for any " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 1 " " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"c G" }{TEXT -1 382 " is the largest root, we conclude that the largest e igenvalue always corresponds to m = 1, is real, and is the unique eige nvalue in the range [0.47+0.2c, 0.52+0.2c].\n\n On steps 1 and 3 we ha ve to show that roots of a polynomial are inside a circle of radius r in the complex plane. This task is similar to the task of establishi ng\nstability of a filter with the transfer function " }{XPPEDIT 18 0 "1/a(z)" "6#*&\"\"\"\"\"\"-%\"aG6#%\"zG!\"\"" }{TEXT -1 8 ", where " } {XPPEDIT 18 0 "a(z)" "6#-%\"aG6#%\"zG" }{TEXT -1 706 " is a polynomial . Such filter is stable, if all roots of the polynomial are inside th e unit circle. \nA variety of tests exist for this condition; for our \+ purposes, the algebraic Marden-Jury test is convenient. With aproporia te rescaling of the variable it can be used to prove that all roots of a polynomial are inside the circle of any given radius r. As the test requires only a simple algebraic calculaion on the coefficients of th e polynomial, it can be easily performed for symbolic and interval coe fficients.\nFinally, we compute the largest root of the characteristic polynomial numerically for all valences up to some maximum. For each computed root, we verify that that the precision is at least " } {XPPEDIT 18 0 "epsilon = 1e-11" "6#/%(epsilonG$\"\"\"!#6" }{TEXT -1 61 ": we use interval arithmetics to evaluate the polynomial at " } {XPPEDIT 18 0 "lambda[0] - epsilon" "6#,&&%'lambdaG6#\"\"!\"\"\"%(epsi lonG!\"\"" }{TEXT -1 7 " and " }{XPPEDIT 18 0 "lambda[0]+epsilon" "6 #,&&%'lambdaG6#\"\"!\"\"\"%(epsilonGF(" }{TEXT -1 270 " and assert tha t the sign is guaranteed to change. There may be more than one root : we still have to prove that there is only a single root in the compu ted interval and that the rest of the roots are smaller. The maximal \+ valence N is chosen in such a way that for N " }{XPPEDIT 18 0 "cos(2*P i/N)" "6#-%$cosG6#*(\"\"#\"\"\"%#PiGF(%\"NG!\"\"" }{TEXT -1 121 "is \" sufficiently close\" to 1. This means that for all K > N correspo nding eigenvalue differs from the limit value " }{XPPEDIT 18 0 "lamb da[infinity]" "6#&%'lambdaG6#%)infinityG" }{TEXT -1 17 " by no more th an " }{XPPEDIT 18 0 "epsilon" "6#%(epsilonG" }{TEXT -1 10 " , where \+ " }{XPPEDIT 18 0 "epsilon" "6#%(epsilonG" }{TEXT -1 190 " is small en ough for us to establish, using interval arithmetics, that the Jacob ian of the characteristic map is positive for eigenvectors computed u sing formulas derived below for all " }{XPPEDIT 18 0 "lambda " "6#%'l ambdaG" }{TEXT -1 17 " in the interval " }{XPPEDIT 18 0 "[lambda[infin ity]-epsilon,lambda[infinity]]" "6#7$,&&%'lambdaG6#%)infinityG\"\"\"%( epsilonG!\"\"&F&6#F(" }{TEXT -1 171 " . The actual computation of the \+ Jacobian and evaluation of the necessary contraction functions is perf ormed in the C part of the code (This is not done yet; the value of " }{XPPEDIT 18 0 "epsilon;" "6#%(epsilonG" }{TEXT -1 78 " can be detrmin ed only after the C code is made to to work with dual schemes)." }} {SECT 0 {PARA 4 "" 0 "" {TEXT -1 16 "Marden-Jury test" }}{PARA 0 "" 0 "" {MPLTEXT 0 21 26 "MardenJury(a,var,rootrad) " }{TEXT -1 44 "Compute Marden-Jury table for a polynomial " }{MPLTEXT 0 21 1 "p" }{TEXT -1 13 " in variable " }{MPLTEXT 0 21 3 "var" }{TEXT -1 117 ", with variab le rescaled by rootrad.\nUsed to verify that all roots of a polynomial are inside the circle of radius r." }}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "MardenJury := proc( a::polynom, var::name, rootrad) loc al i,k,M, acol, tbl, restable;\n M := degree(a,var);\n for i from \+ 0 to M do \n tbl[0,i] := coeff(a, var,i)*rootrad^(i-M);\n od;\n \+ for i from 1 to M do\n for k from 0 to M-i do\n tbl[i,k] := tbl[i-1,0]*tbl[i-1,k]- tbl[i-1,M-k-i+1]*tbl[i-1,M-i+1];\n od\n \+ od;\n for i from 1 to M do \n restable[i] := tbl[i,0]; \n od; \n eval(restable);\nend: " "6#>%+MardenJuryGR6%'%\"aG%(polyno mG'%$varG%%nameG%(rootradG7(%\"iG%\"kG%\"MG%%acolG%$tblG%)restableG6\" F5C'>F1-%'degreeG6$F(F+?(F/\"\"!\"\"\"F1%%trueG>&F36$F?(F0F>&F36$F/ F0,&*&&F36$,&F/F=\"\"\"FHF>&F46#F/&F36$F/F<-%%evalG6#F4F5F5F5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 36 "Interval version of Marden-Jury test" }}{PARA 0 "> " 0 "" {XPPEDIT 19 1 "IntervMardenJury := proc( a::polynom(interv al), var::name, rootrad::numeric) local i,k,M, acol, tbl, restable;\n \+ # checking correct type of a is too messy\n M := degree(a,var);\n \+ for i from 0 to M do \n tbl[0,i] := Interval_times(coeff(a, var, i), rootrad^(i-M));\n od;\n for i from 1 to M do\n for k from \+ 0 to M-i do\n tbl[i,k] := Interval_add( Interval_times( tbl[i-1, 0], tbl[i-1,k]), \n Interval_times(-1, Interval_times( tbl[i- 1,M-k-i+1], tbl[i-1,M-i+1]) ));\n od\n od;\n for i from 1 to M \+ do \n restable[i] := tbl[i,0]; \n od;\n eval(restable);\nend: " "6#>%1IntervMardenJuryGR6%'%\"aG-%(polynomG6#%)intervalG'%$v arG%%nameG'%(rootradG%(numericG7(%\"iG%\"kG%\"MG%%acolG%$tblG%)restabl eG6\"F:C'>F6-%'degreeG6$F(F.?(F4\"\"!\"\"\"F6%%trueG>&F86$FAF4-%/Inter val_timesG6$-%&coeffG6%F(F.F4)F1,&F4FBF6!\"\"?(F4\"\"\"FBF6FC?(F5FAFB, &F6FBF4FOFC>&F86$F4F5-%-Interval_addG6$-FH6$&F86$,&F4FB\"\"\"FOFA&F86$ ,&F4FB\"\"\"FOF5-FH6$,$\"\"\"FO-FH6$&F86$,&F4FB\"\"\"FO,*F6FBF5FOF4FO \"\"\"FB&F86$,&F4FB\"\"\"FO,(F6FBF4FO\"\"\"FB?(F4\"\"\"FBF6FC>&F96#F4& F86$F4FA-%%evalG6#F9F:F:F:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}} {SECT 0 {PARA 4 "" 0 "" {TEXT -1 9 "Deflation" }}{PARA 0 "" 0 "" {MPLTEXT 0 21 23 "deflate(p,var,rootval) " }{TEXT -1 43 "compute the c oefficients of the polynomial " }{XPPEDIT 18 0 "p(z)/(z -z[0])" "6#*&- %\"pG6#%\"zG\"\"\",&F'F(&F'6#\"\"!!\"\"F-" }{TEXT -1 39 "; it is assum ed that p is divisible by " }{XPPEDIT 18 0 "z-z[0]" "6#,&%\"zG\"\"\"&F $6#\"\"!!\"\"" }{TEXT -1 3 ". " }{MPLTEXT 0 21 3 "var" }{TEXT -1 30 " is the name of the variable, " }{MPLTEXT 0 21 7 "rootval" }{TEXT -1 13 " is the root." }}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "deflate : = proc(p::polynom, var::name, rootval)\nlocal i, dp,r; \n dp := 0; \n r := lcoeff(p,var);\n for i from degree(p,var)-1 to 0 by -1 d o \n dp := dp + r*var^i;\n r := coeff(p, var, i ) + rootval *r;\n od; \n dp;\nend: \n" "6#>%(deflateGR6%'%\"pG%(polynomG'% $varG%%nameG%(rootvalG7%%\"iG%#dpG%\"rG6\"F2C&>F0\"\"!>F1-%'lcoeffG6$F (F+?(F/,&-%'degreeG6$F(F+\"\"\"\"\"\"!\"\",$\"\"\"FAF5%%trueGC$>F0,&F0 F?*&F1F?)F+F/F?F?>F1,&-%&coeffG6%F(F+F/F?*&F-F?F1F?F?F0F2F2F2" }}}} {SECT 0 {PARA 4 "" 0 "" {TEXT -1 27 "Analysis of the eigenvalues" }} {PARA 0 "" 0 "" {TEXT -1 41 "Now we perform steps 1-3 described above. " }}{SECT 0 {PARA 0 "" 0 "" {TEXT -1 73 "(1). We show that for c < 0, \+ all roots of the characteristic polynomial " }{XPPEDIT 18 0 "P(c,lamb da)" "6#-%\"PG6$%\"cG%'lambdaG" }{TEXT -1 142 " are less than 0.51 (ac tually, they are less than 0.5, but due to numerical nature of our cal culations, we have to relax the upper boundary). " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "MJtab := MardenJury(CharPolyMainFactor, lambda, 51/100):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "MJtabInterv := map( unapply( 'inapp ly'( dummy, c), dummy), MJtab ):" }}}{EXCHG {PARA 0 "> " 0 "" {XPPEDIT 19 1 "TestNegativeC := proc( cstart::numeric, cend::numeric, cstep::numeric ) \n local MJ,cx; global MJtabInterv;\n for cx from c start to cend by cstep do \n MJ := map( unapply( dummy([cx,cx+cstep ]),dummy), MJtabInterv);\n if op(2, MJ[1]) > 0 or op(1, MJ[2]) < \+ 0 or op(1,MJ[3]) < 0 or op(1, MJ[4]) < 0 \n then ERROR(`test fa iled for interval = `, [cx,cx+cstep]); \n fi;\n od;\n print(`All \+ tests passed`);\nend:" "6#>%.TestNegativeCGR6%'%'cstartG%(numericG'%%c endGF)'%&cstepGF)7$%#MJG%#cxG6\"F1C$?(F0F(F-F+%%trueGC$>F/-%$mapG6$-%( unapplyG6$-%&dummyG6#7$F0,&F0\"\"\"F-FBF>%,MJtabIntervG@$5552\"\"!-%#o pG6$\"\"#&F/6#\"\"\"2-FK6$\"\"\"&F/6#\"\"#FI2-FK6$\"\"\"&F/6#\"\"$FI2- FK6$\"\"\"&F/6#\"\"%FI-%&ERRORG6$% " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 148 "Do tests, ad justing the step in c. This is not really necessary -- we can simply t ake the smallest step, but to save time we use larger steps first." } {MPLTEXT 1 0 32 "\nTestNegativeC(-1.0,-0.65,0.05);" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#%1All~tests~passedG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "TestNegativeC(-0.6,-0.275,0.025);" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#%1All~tests~passedG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "TestNegativeC(-0.25,-0.06,0.01);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%1All~tests~passedG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "TestNegativeC(-0.05,0.,0.005);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%1All~tests~passedG" }}}}{SECT 0 {PARA 0 "" 0 "" {TEXT -1 26 "(2). We show that for any " }{XPPEDIT 18 0 "c = 0..1" "6#/%\"cG ;\"\"!\"\"\"" }{TEXT -1 30 ", there is a unique real root " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 59 " in the interval [0.47+0.2*c,0.52+0. 2*c], and the function " }{XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" } {TEXT -1 12 " increases." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 11 "Deriv ative " }{XPPEDIT 18 0 "diff(lambda,c);" "6#-%%diffG6$%'lambdaG%\"cG" }{TEXT -1 21 " can be computed as " }{XPPEDIT 18 0 "-diff(F,c)/diff(F ,lambda);" "6#,$*&-%%diffG6$%\"FG%\"cG\"\"\"-F&6$F(%'lambdaG!\"\"F." } {TEXT -1 5 " if " }{XPPEDIT 18 0 "F(lambda,c);" "6#-%\"FG6$%'lambdaG% \"cG" }{TEXT -1 67 " is the char. polynomial. We use the derivative o nly to show that " }{XPPEDIT 18 0 "lambda(c);" "6#-%'lambdaG6#%\"cG" } {TEXT -1 159 " increases, thus, we need only the sign. We will see tha t in the domain of interest the numerator is always negative, thus we \+ need to look only on the sign of " }{XPPEDIT 18 0 "diff(F,lambda);" "6 #-%%diffG6$%\"FG%'lambdaG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "Flambda := diff( CharPolyMainFactor, lambda );" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#>%(FlambdaG,0*$)%'lambdaG\"\"%\"\"\"\"\"&*&,&#!#D\"#K \"\"\"%\"cG#!\"&F0F1)F(\"\"$F*F)*&,(F2#\"#<\"$7&#\"#@\"$G\"F1*$)F2\"\" #F*#!\"\"F;F1)F(FAF*F6*&,(F?#\"\"(\"%#>)#F/\"%[?F1F2#!#8FIF1F(F1FAF2#F C\"&oF$#F1\"%'4%F1F?FN" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "F c := diff( CharPolyMainFactor, c);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# >%#FcG,0*$)%'lambdaG\"\"%\"\"\"#!\"&\"#K*&,&#\"#<\"$7&\"\"\"%\"cG#!\" \"\"$c#F3)F(\"\"$F*F3*&,&F4#\"\"(\"%'4%#!#8\"%#>)F3F3)F(\"\"#F*F3*&,&# F6\"&oF$F3F4#F6\"&%Q;F3F(F3F3#F3\"')GC&F3F4#\"\"&\"(w&[5*$)F4FCF*#F9\" (_r4#" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 218 "This plot illustrates the idea: for c > 0, the largest root can be bounded from above and below by linear functions of c (green lines). The area between the lines is above the red cur ve \nindicating the zero set of " }{XPPEDIT 18 0 "diff(F,lambda);" "6 #-%%diffG6$%\"FG%'lambdaG" }{TEXT -1 31 ", where derivative is positiv e." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 240 "display( contourplot( unappl y( Flambda,[c, lambda]), -1..1, 0..1, contours=[0], color=red), seq(pl ot(op(i,[EigenvalsPlotLists]),color=blue), i = 1..5), plot( 0.47+0.2*c , c = -1..1, color=green), plot( 0.52+0.2*c, c = -1..1, color = green) );" }}{PARA 13 "" 1 "" {GLPLOT2D 701 480 480 {PLOTDATA 2 "6,-%'CURVESG 6cx7$7$$!\"\"\"\"!$\"1Uz?[JD`J!#<7$$!1*[c&*)*p=$)*!#;$\"1VCyZ*\\$fJF-7 $7$$!1-++++++#*F1$\"1^j?.c%f3$F-F.7$7$F($\"1#eG9dGRy&F-7$$!1sO=HH3y&*F 1$\"1d$=fk9/*eF-7$7$$!1+++++++#*F1$\"1\\z1&=GC*fF-F>7$7$FE$\"19v'4B>[f \"F17$$!1BQ1x'[-B*F1$\"1+++++++;F17$FM7$F($\"1F&*o*F1$\"1h5 \\*f:$4EF1FZ7$F57$$!17qp)3X#=!*F1$\"1b][VaA\"4$F-7$7$$!1,++++++%)F1$\" 1_]]rYyCIF-F`o7$FD7$$!1Z-g!HLp\"))F1$\"1N7+`km%3'F-7$7$Fgo$\"1mH_Q2j%= 'F-F\\p7$FJ7$$!1&G#45cY)=*F1$\"1Uh/0GB%f\"F17$7$Fgo$\"1qfWET-([\"F1Ffp 7$7$FE$\"1i5\\*f:$4EF17$$!1)[U@W>!>*)F1$\"1W72@(4&fEF17$7$$!1A)4ChG#4' )F1$\"1+++++++GF1Fcq7$Fiq7$$!11y-_I+M%)F1$\"1.R,E:+!> 8Fhx9F17$7 $F`s$\"1#oCCajaT\"F1Fdt7$Fer7$$!1'*R9?(QJy(F1$\"1)*>2g$p:*GF17$7$F_t$ \"12*e@$\\ZIHF1F^u7$F_s7$$!10?G4/V%R(F1$\"1:+TY?:sHF-7$7$$!1-++++++oF1 $\"1!e*y4/Y[C&G(F1$\"1D%3&)4CiU 'F-7$7$F_v$\"1=OsSY]7lF-Fgv7$Fjt7$$!1PK#fngt?(F1$\"1=;'zL!o.9F17$7$F_v $\"1fFr9&z^O\"F1Faw7$7$F`s$\"11*e@$\\ZIHF17$$!1%pB-eA<<(F1$\"1Y=6!Hhe) HF17$7$F_v$\"1V(3*4Su!3$F1F^x7$7$F_v$\"1z&*y4/YHF-7$7$$!1-++++++gF1$\"1Gd9kx$*pGF-F[y7$F]w7$$!1&)*GKh6S^'F1$ \"1:\\9m!e+d'F-7$7$Fby$\"1:&Q0i(RXmF-Fgy7$Fgw7$$!1:LB!G\"z1jF1$\"1cm6S cR`8F17$7$Fby$\"16ajl@XG8F1Faz7$Fdx7$$!1())3pA(*zh'F1$\"1UWX8')**3JF17 $7$$!1n(*eMuMLjF1$\"1+++++++KF1F[[l7$Fa[l7$$!1V&o71lP0'F1$\"1rUjID)oA$ F17$7$Fby$\"1B&y4c8dB$F1Fg[l7$Fay7$$!1.Ez>jv6D#f+I$F17$7$Fh \\l$\"1(yu;7SsL$F1Fa^l7$Fg\\l7$$!1`Q?r@4l\\F1$\"1_#>g&3YDGF-7$7$$!1-++ ++++WF1$\"15Q**pRb%y#F-F[_l7$Fc]l7$$!1OzCu/ff\\F1$\"1m'R7P_zz'F-7$7$$! 1.++++++WF1$\"1$=P?:6c%oF-Fg_l7$F]^l7$$!15EQ3b]#e%F1$\"1/8>aFD\"H\"F17 $7$Fb_l$\"1E$[&Gt#)z7F1Fc`l7$Fg^l7$$!1%p>@bmqx%F1$\"1Y)fgFL&)Q$F17$7$F ^`l$\"1Q9U\"*\\:pMF1F]al7$Fa_l7$$!1d)H^,#ecTF1$\"1w#\\c25Hy#F-7$7$$!1- ++++++OF1$\"1CxgwI(eu#F-Fgal7$7$Fb_lF``l7$$!1uzUC1=wTF1$\"1h)R@7.4)oF- 7$7$F^bl$\"10YJ[8#H\"pF-Fdbl7$Fi`l7$$!1hh&**>!*Hu$F1$\"1z!y**4&\\r7F17 $7$F^bl$\"1gId=Z_j7F1F^cl7$Fcal7$$!19RL5wR&>%F1$\"1cp;0))p(\\$F17$7$$! 1?*>jHH!*y$F1$\"1+++++++OF1Fhcl7$F^dl7$$!1;+m!*G=IOF1$\"11+LX94:OF17$7 $F^bl$\"1%y()))\\b'>OF1Fddl7$F]bl7$$!1&z+clq&[LF1$\"1lR+yK&Gu#F-7$7$$! 1-++++++GF1$\"1pS&fI7%4FF-F^el7$Fjbl7$$!1J=t.#H')Q$F1$\"1U\"f'=g9VpF-7 $7$Feel$\"1=-%p\\f!fpF-Fjel7$Fdcl7$$!1$G\"oIv$H\"HF1$\"1S1Ml(okD\"F17$ 7$Feel$\"1Hfk?B#3D\"F1Fdfl7$Fjdl7$$!1s6rIiWqHF1$\"1&eb`6B_o$F17$7$Feel $\"1'fl'4Pu8PF1F^gl7$Fdel7$$!1^T\\+E*4a#F1$\"1Y2Z-I'\\q#F-7$7$$!1-++++ ++?F1$\"1KBRK])[n#F-Fhgl7$F`fl7$$!1$y>)))o-(f#F1$\"14*)4WW8&)pF-7$7$F_ hl$\"1uT?[0*\\)pF-Fdhl7$Fjfl7$$!1%Qy>()G**3#F1$\"1\">*)fVk\\C\"F17$7$F _hl$\"1E[SXr!4C\"F1F^il7$Fdgl7$$!1FXp\\SrMBF1$\"1js%[-dtw$F17$7$F_hl$ \"1$yRE%[.IQF1Fhil7$F^hl7$$!1/r[0-zL`,=F1$\"1+g$phfw+(F-7$7$Fijl$\"1]' p`Ve>*pF-F^[m7$Fdil7$$!1'=h2'eIs7F1$\"1#f!QIH:O7F17$7$Fijl$\"1L@(G:%>L 7F1Fh[m7$F^jl7$$!1(fM&QGvHk(['QF17$7$Fijl$\"1%o8rJnu(RF1Fb \\m7$Fhjl7$$!1x3lY;7p#*F-$\"1HaKB3cMEF-7$7$$!1<++++++SF-$\"1ubPzZk5EF- F\\]m7$7$Fijl$\"1['p`Ve>*pF-7$$!1-z$y/fB+\"F1$\"1,&*=R_z6qF-7$7$Fc]m$ \"1CrqfkT\")pF-F[^m7$F^\\m7$$!16Z>1W$*)e%F-$\"1N(4.sY%H7F17$7$Fc]m$\"1 %RVXf^sA\"F1Fe^m7$Fh\\m7$$!1#*[*y$e8l6F1$\"1Xu%*=zc#)RF17$7$$!1)\\nJwr T4\"F1$\"1+++++++SF1F__m7$7$Ff_m$\"1,++++++SF17$$!1wI_\"*)RC<&F-$\"1`h d%*>ieSF17$7$Fc]m$\"1Kkjc4pwSF1F^`m7$Fb]m7$$!1%=UniAJ?\"F-$\"1%3rLJh:g #F-7$7$$\"1$)************RF-$\"1Er0N3]!e#F-Fh`m7$Fa^m7$$!13#RV%Hn(*>F- $\"1&fp@ZO))*pF-7$7$F_am$\"1(4L!4*e\\&pF-Fdam7$F[_m7$$\"1ZLzjCk6NF-$\" 1K.\"o(yTC7F17$7$F_am$\"17P#R*QvA7F1F^bm7$Fd`m7$$\"1SyebHF%H\"F-$\"121 A_jGNTF17$7$F_am$\"1'*3`Eq1\"=%F1Fhbm7$F^am7$$\"1zW)ow;0'oF-$\"1^xb;;u pDF-7$7$$\"1)************>\"F1$\"1WL^j'=9b#F-Fbcm7$Fjam7$$\"1&=z7pC&fg F-$\"1*RgVlP-(pF-7$7$Ficm$\"1)R\\![wD9pF-F^dm7$Fdbm7$$\"1W(f?N0&e6F1$ \"1F,(RKZ2A\"F17$7$Ficm$\"1t8,,4X>7F1Fhdm7$F^cm7$$\"1B&pv?\\E_(F-$\"1A :iRv'QA%F17$7$Ficm$\"1F],%)z#yI%F1Fbem7$Fhcm7$$\"1am['p=A\\\"F1$\"1@nc F1$\"1@\"fk1CK_#F-F\\fm7$Fddm7$$\"1C^O' \\*[99F1$\"1rV<=DbFpF-7$7$Fcfm$\"1(=3$*G15'oF-Fhfm7$F^em7$$\"1^.A0ahj> F1$\"1C)*Q(H#>=7F17$7$Fcfm$\"1@-!fCZr@\"F1Fbgm7$Fhem7$$\"1ZPUWJBX8F1$ \"1D\")yFMQFVF17$7$$\"16m`O/h!p\"F1$\"1*************R%F1F\\hm7$Fbhm7$$ \"1Tk1)yC'[>F1$\"1xn'fg(oDWF17$7$Fcfm$\"1pv15;5LWF1Fhhm7$7$Fcfm$\"1?\" fk1CK_#F-7$$\"1n4b3bA)H#F1$\"1e^CdC()3DF-7$7$$\"1*************z#F1$\"1 DM+AEv&\\#F-Feim7$F^gm7$$\"137HWS_DAF1$\"1^Ray(zB(oF-7$7$F\\jm$\"1?$Ga Wtoz'F-Fajm7$7$Fcfm$\"1?-!fCZr@\"F17$$\"1?Z2k.(ow#F1$\"1SE'z\"[c;7F17$ 7$F\\jm$\"1WF>A\")o:7F1F^[n7$F^im7$$\"1w5S+sD/EF1$\"1g%*z*Rry\\%F17$7$ F\\jm$\"1e#*\\u&f%GXF1Fh[n7$7$F\\jm$\"1CM+AEv&\\#F-7$$\"1&\\*3Es5/JF1$ \"1@DbpQYzCF-7$7$$\"1*************f$F1$\"17>ogk%)oCF-Fe\\n7$Fgjm7$$\"1 f%p#)zO(QIF1$\"1.Fl3gJ1oF-7$7$F\\]n$\"1O[t$flMs'F-Fa]n7$Fd[n7$$\"12upt ecoNF1$\"1'H^J1>='zXF17$7$F\\]n$\"1bg-yd\"3k%F1Fe^n7$F[]n7$$\"1SW%z0)*) 4RF1$\"1$zx-r40X#F-7$7$Fehm$\"1lOsv6NUCF-F__n7$7$F\\]n$\"1P[t$flMs'F-7 $$\"1nN)zj;Q&QF1$\"1f@35o\"4t'F-7$7$Fehm$\"1axD>$*HUmF-F\\`n7$7$$\"1)* ***********f$F1Fb^n7$$\"1Pn(RiP*oVF1$\"1J;,)=Jb@\"F17$7$Fehm$\"1OkgQ)H [@\"F1Fi`n7$F[_n7$$\"1npWHG2aQF1$\"1:lF&ejHn%F17$7$Fehm$\"1-7,41;vZF1F can7$7$$\"1)************R%F1Ff_n7$$\"1%)*GQ*>j:ZF1$\"1!3b3.S=U#F-7$7$$ \"1*************>&F1$\"1IJMPF6;CF-F`bn7$Fb`n7$$\"1\\#yzpm/n%F1$\"1X(3, ^mwk'F-7$7$$\"1+++++++_F1$\"1#4k5;#yalF-F\\cn7$F_an7$$\"1;\"**Q\"eT$\\F17$7$Fbfn$\" 1$>[#\\%yk(\\F1F^hn7$Fafn7$$\"1IJBJ$oqK'F1$\"1WV$QMeYO#F-7$7$$\"1)**** ********z'F1$\"10-iwDxjBF-Fhhn7$7$Fbfn$\"1+y$p&4?ikF-7$$\"1k?)GQ'Q2jF1 $\"1q'*e&3oIY'F-7$7$F_in$\"1us5HkrljF-Fgin7$Fjgn7$$\"1>8V2\\%Qw'F1$\"1 RVGYv2=7F17$7$F_in$\"1G1&yeGu@\"F1Fajn7$Fdhn7$$\"17[iQH=gjF1$\"1\"f(oI &3*>]F17$7$$\"1(************z'F1$\"1/#3h;]^4&F1F[[o7$F^in7$$\"1'z%GyU% G8(F1$\"16gd3'ydL#F-7$7$$\"1)************f(F1$\"1p6))=TLPBF-Fg[o7$F]jn 7$$\"19TOAi;FrF1$\"1A%z\")))oTO'F-7$7$F^\\o$\"1(pGrtijE'F-Fc\\o7$Fgjn7 $$\"1/+Hv*40c(F1$\"1'*\\N7]u>7F17$7$F^\\o$\"1Bs,Fv3>7F1F]]o7$Fa[o7$$\" 1`@#p(=EmpF1$\"1A*Q:1po6&F17$7$$\"1T+W6-x6uF1FgbnFg]o7$F]^o7$$\"1\\(>H gR!ovF1$\"1A,a)>!)f@&F17$7$F^\\o$\"1EXA'eF/A&F1Fa^o7$F]\\o7$$\"1f;%RE7 (QzF1$\"1$p\"H!oQkI#F-7$7$$\"1(************R)F1$\"1L!=eSq/J#F-F[_o7$Fi \\o7$$\"1C\"pi,Rv%zF1$\"1oVl=\\IiiF-7$7$$\"1'************R)F1$\"17<7,3 0lhF-Fg_o7$7$$\"1(************f(F1Fd]o7$$\"10SwF*>lN)F1$\"1'*z6O+u@7F1 7$7$Fb_o$\"1pwgc$o5A\"F1Ff`o7$7$Fd`o$\"1DXA'eF/A&F17$$\"1m9)*oabB#)F1$ \"1l#4bEA#)G&F17$7$Fb_o$\"1JRuSnU9`F1Fcao7$Fa_o7$$\"1&o5T[*F1$\"1l:13B(HG#F-F]bo7$7$Fb_oF``o7$$\"1 E#p^?>$o()F1$\"1dQ:uRSehF-7$7$Fdbo$\"1`-OKVdigF-Fjbo7$7$F^`oF]ao7$$\"1 <\"[8,]>:*F1$\"1TfK%*\\-C7F17$7$Fdbo$\"1#>zeBOLA\"F1Feco7$Fiao7$$\"1Q% HzDiN'))F1$\"1z_.r)=#o`F17$7$Fdbo$\"1pU,31y@aF1F_do7$Fcbo7$$\"1!f5k*3 \"4b*F1$\"1Qq%z^XaC#F-7$7$$\"\"\"F*$\"1`'zu*HgaAF-Fido7$F`co7$$\"1>M)e 9W$*e*F1$\"1'*Geq#zK0'F-7$7$F`eo$\"1WZHpwhffF-Feeo7$F[do7$$\"1!G)Q3U'o %**F1$\"1ge!e*ycE7F17$7$F`eo$\"1_;Wb7'eA\"F1F_fo7$Fedo7$$\"1!))pT_Ea[* F1$\"1f]\"zt'GdaF17$7$F`eo$\"1K%yp&yaXbF1Fifo-%'COLOURG6&%$RGBG$\"*+++ +\"!\")F*F*-F$6$7`q7$F($\"1++++++#>)!#G7$$!1+++++++)*F1$\"1+++RFle6!#> 7$$!1+++++++'*F1$\"1+++TL?!H#Feho7$$!1*************R*F1$\"1+++t=j$R$Fe ho7$FE$\"1+++8O(yY%Feho7$$!1+++++++!*F1$\"1+++c$>=^&Feho7$$!1+++++++)) F1$\"1+++I^JClFeho7$$!1+++++++')F1$\"1+++:J;/vFeho7$$!1+++++++%)F1$\"1 +++V?7]%)Feho7$$!1+++++++#)F1$\"1+++#>34O*Feho7$$!1+++++++!)F1$\"1++++ (>N-\"!#=7$$!1+++++++yF1$\"1+++[C;26Fa[p7$$!1+++++++wF1$\"1+++M*yo=\"F a[p7$$!1+++++++uF1$\"1+++\"eDDE\"Fa[p7$$!1+++++++sF1$\"1+++Lc&RL\"Fa[p 7$$!1+++++++qF1$\"1+++;&>5S\"Fa[p7$$!1+++++++oF1$\"1+++e_cj9Fa[p7$$!1+ ++++++mF1$\"1+++H*Q9_\"Fa[p7$$!1+++++++kF1$\"1+++e^[u:Fa[p7$$!1+++++++ iF1$\"1+++:xaA;Fa[p7$$!1+++++++gF1$\"1+++W,Zl;Fa[p7$$!1+++++++eF1$\"1+ ++Ek4.s9Fa[p7$$!1+++++++CF1$\"1+++\"zK/S \"Fa[p7$$!1+++++++AF1$\"1+++'34;K\"Fa[p7$$!1+++++++?F1$\"1+++<)GdB\"Fa [p7$$!1+++++++=F1$\"1+++cO\"G9\"Fa[p7$$!1+++++++;F1$\"1+++Q*)*G/\"Fa[p 7$$!1+++++++9F1$\"1+++#[L.O*Feho7$$!1+++++++7F1$\"1+++!*\\yA#)Feho7$$! 1+++++++5F1$\"1******yR3YKFeho7 $$\"1+++++++gF-$\"1+++:(G!e\\Feho7$$\"1+++++++!)F-$\"1+++^!\\qs'Feho7$ $\"1+++++++5F1$\"1+++7\"=pZ JtFa[p7$$\"1+++++++mF1$\"1+++^p#4h(Fa[p7$$\"1+++++++oF1$\"1,++[%*>#*yF a[p7$$\"1+++++++qF1$\"1+++s%=_<)Fa[p7$$\"1+++++++sF1$\"1*****zI6*f%)Fa [p7$$\"1+++++++uF1$\"1+++\"R3iu)Fa[p7$$\"1+++++++wF1$\"1+++tJ/M!*Fa[p7 $$\"1+++++++yF1$\"1+++4?NB$*Fa[p7$$FehpF1$\"1+++`S29'*Fa[p7$$\"1++++++ +#)F1$\"1+++U6:1**Fa[p7$$\"1+++++++%)F1$\"1+++oF&*>5F-7$$\"1+++++++')F 1$\"1+++b]T\\5F-7$$\"1+++++++))F1$\"1+++)*o**y5F-7$$\"1+++++++!*F1$\"1 +++[Mp36F-7$$\"1+++++++#*F1$\"1+++r+]Q6F-7$$\"1*************R*F1$\"1++ +SBTo6F-7$$\"1+++++++'*F1$\"1+++DgU)>\"F-7$$\"1+++++++)*F1$\"1+++(3P&G 7F--Fcgo6&FegoF*F*Ffgo-F$6$7`q7$F($\"1+++(*****\\iF-7$Faho$\"1+++/=tfb F-7$Fgho$\"1+++qzG#H&F-7$F\\io$\"1+++)HIW4&F-7$FE$\"1+++'y8E$\\F-7$Fdi o$\"1+++*)p&Qz%F-7$Fiio$\"1+++LD^rYF-7$F^jo$\"1+++U0khXF-7$Fcjo$\"1+++ r8%F-7$F]\\p$\"1+++U*3E1%F-7$Fb\\p$\"1+++C=P(*RF-7$ Fg\\p$\"1+++d.&e$RF-7$F\\]p$\"1+++X$Gx(QF-7$Fa]p$\"1+++=@uAQF-7$Ff]p$ \"1+++8)p1x$F-7$F[^p$\"1+++uLK@PF-7$F`^p$\"1+++VAauOF-7$Fe^p$\"1+++A') =IOF-7$Fj^p$\"1+++7O9)e$F-7$F__p$\"1+++eTI[NF-7$Fd_p$\"1+++Q0e5NF-7$Fi _p$\"1+++BV*[Z$F-7$F^`p$\"1+++xmJp,KF-7$Fecp$\"1+++[Jv'=$F-7$Fjcp$\"1+++pl MtJF-7$F_dp$\"1+++_AXhJF-7$Fddp$\"1+++$*)\\5:$F-7$Fidp$\"1+++j)>@9$F-7 $F^ep$\"1+++*3VY8$F-7$Fcep$\"1+++15gGJF-7$Fhep$\"1+++Ia(R7$F-7$F]fp$\" 1+++z&[27$F-7$Fbfp$\"1+++XI!*=JF-7$Fffp$\"1+++0=U=JF-7$Fjfp$\"1+++V#)G >JF-7$F^gp$\"1+++%>'[@JF-7$F*$\"1++++++DJF-7$Ffgp$\"1+++5Y\")HJF-7$F[h p$\"1+++3c\"f8$F-7$F_hp$\"1+++R$*GVJF-7$Fdhp$\"1+++,I#>:$F-7$Fihp$\"1+ ++VZ!=;$F-7$F^ip$\"1+++pP#H<$F-7$Fcip$\"1+++b/F&=$F-7$FP$\"1+++hk$))>$ F-7$F[jp$\"1+++!*[h8KF-7$F`jp$\"1+++m.gHKF-7$Fdjp$\"1+++.#*yYKF-7$Fijp $\"1+++]&z^E$F-7$F^[q$\"1+++Q:x%G$F-7$F\\r$\"1+++Nuc0LF-7$Ff[q$\"1+++M =dFLF-7$Fd[l$\"1+++K=z]LF-7$F^\\q$\"1+++)GP_P$F-7$Fadl$\"1+++'*4#4S$F- 7$Ff\\q$\"1+++***eyU$F-7$Fh_m$\"1+++x32cMF-7$F^]q$\"1+++'4!e&[$F-7$Fc] q$\"1+++%R9k^$F-7$Fh]q$\"1+++Qig[NF-7$F]^q$\"1+++fL>#e$F-7$Fb^q$\"1+++ ;%>sh$F-7$Fccn$\"1+++sYt`OF-7$Fj^q$\"1+++Ipz\"p$F-7$F__q$\"1+++uDZJPF- 7$Fd_q$\"1+++my$Gx$F-7$Fi_q$\"1+++y/)f\"QF-7$F]`q$\"1+++O9+hQF-7$Fb`q$ \"1+++*[TF-7$F[bq$\"1+++r!p'yTF-7$F`bq$\"1++++q\\TUF-7 $Febq$\"1+++*H!)zI%F-7$Fjbq$\"1+++IifyVF-7$F^cq$\"1+++hs%RX%F-7$Fccq$ \"1+++,3\"[`%F-7$Fhcq$\"1+++y!=Ai%F-7$F]dq$\"1+++r`ey]F-7$Faeq$\"1+++%yXKC& F-7$Ffeq$\"1+++yCXhaF-Fjeq-F$6$7`qF_fq7$Faho$\"1+++8(eD+(F-7$Fgho$\"1, ++([!*=L(F-7$F\\io$\"1+++2zD\"f(F-7$FE$\"1*****\\$4B9yF-7$Fdio$\"1**** *>/;Q,)F-7$Fiio$\"1,++*3'o'>)F-7$F^jo$\"1+++\"[1oO)F-7$Fcjo$\"1+++MDwE &)F-7$Fhjo$\"1+++t5Ny')F-7$F][p$\"1+++O;)G#))F-7$Fc[p$\"1+++KMMh*)F-7$ Fh[p$\"1+++wi]%4*F-7$F]\\p$\"1+++\")G)HA*F-7$Fb\\p$\"1+++65FZ$*F-7$Fg \\p$\"1+++&H\"yn%*F-7$F\\]p$\"1+++jy&[e*F-7$Fa]p$\"1+++HAz)p*F-7$Ff]p$ \"1+++)RM)4)*F-7$F[^p$\"1+++$>+#=**F-7$F`^p$\"1+++OyS-5F17$Fe^p$\"1+++ RMw75F17$Fj^p$\"1+++0:!H-\"F17$F__p$\"1+++e]$G.\"F17$Fd_p$\"1+++.edU5F 17$Fi_p$\"1+++XU8_5F17$F^`p$\"1+++;*>:1\"F17$Fc`p$\"1+++>9uq5F17$Fh`p$ \"1+++hm!)z5F17$F]ap$\"1+++DGs)3\"F17$Fbap$\"1+++^l\\(4\"F17$Fgap$\"1+ +++R816F17$F\\bp$\"1+++J0k96F17$Fabp$\"1+++X;-B6F17$Ffbp$\"1+++D@GJ6F1 7$F[cp$\"1+++SkUR6F17$F`cp$\"1+++^)eu9\"F17$Fecp$\"1+++sKQb6F17$Fjcp$ \"1+++NM?j6F17$F_dp$\"1+++eF#4<\"F17$Fddp$\"1+++[Xay6F17$Fidp$\"1+++c= 2'=\"F17$F^ep$\"1+++]w]$>\"F17$Fcep$\"1+++!fa3?\"F17$Fhep$\"1+++v`637F 17$F]fp$\"1+++2CH:7F17$Fbfp$\"1+++E\")QA7F17$Fffp$\"1+++7YSH7F17$Fjfp$ \"1+++,TMO7F17$F^gp$\"1+++6%3KC\"F17$F*$\"1+++0(***\\7F17$Ffgp$\"1+++, ;7R7F17$F[hp$\"1+++D!G$G7F17$F_hp$\"1+++(f8w@\"F17$Fdhp$\"1+++1N(p?\"F 17$Fihp$\"1+++bES'>\"F17$F^ip$\"1+++@j*e=\"F17$Fcip$\"1+++w'\\a<\"F17$ FP$\"1+++m!e];\"F17$F[jp$\"1+++1nra6F17$F`jp$\"1+++**4UW6F17$Fdjp$\"1+ ++Ai;M6F17$Fijp$\"1+++:x%R7\"F17$F^[q$\"1+++'pgP6\"F17$F\\r$\"1+++1/g. 6F17$Ff[q$\"1+++e>Y$4\"F17$Fd[l$\"1+++3/M$3\"F17$F^\\q$\"1+++91Bt5F17$ Fadl$\"1+++ct7j5F17$Ff\\q$\"1+++q^-`5F17$Fh_m$\"1+++9%=H/\"F17$F^]q$\" 1+++V6!G.\"F17$Fc]q$\"1+++armA5F17$Fh]q$\"1+++q)4D,\"F17$F]^q$\"1+++HB K-5F17$Fb^q$\"1+++T0(4#**F-7$Fccn$\"1+++K1E=)*F-7$Fj^q$\"1+++Jp+:(*F-7 $F__q$\"1+++6b66'*F-7$Fd_q$\"1,++*G$[1&*F-7$Fi_q$\"1+++>q*4S*F-7$F]`q$ \"1+++i3`%H*F-7$Fb`q$\"1+++@V%p=*F-7$Fg`q$\"1*****\\sy!y!*F-7$F\\aq$\" 1+++^Svn*)F-7$Faaq$\"1,++yLwb))F-7$Ffaq$\"1+++7p'=u)F-7$F[bq$\"1+++pIy D')F-7$F`bq$\"1+++#Rxr])F-7$Febq$\"1+++(HYcQ)F-7$Fjbq$\"1,++vZpg#)F-7$ F^cq$\"1*****z-.<8)F-7$Fccq$\"1+++^r(y*zF-7$Fhcq$\"1,++<)p\"eyF-7$F]dq $\"1+++D$\\6r(F-7$Fbdq$\"1+++`^wavF-7$Fgdq$\"1+++s;*eQ(F-7$F\\eq$\"1** ***pX=$*>(F-7$Faeq$\"1+++*e_])pF-7$Ffeq$\"1*****zH7#=nF-Fjeq-F$6$7`q7$ F($\"1+++,+++DF17$Faho$\"1+++:%[AG#F17$Fgho$\"1+++EE,)>#F17$F\\io$\"1+ ++k/ZN@F17$FE$\"1+++vr0%3#F17$Fdio$\"1+++@apR?F17$Fiio$\"1+++37I+?F17$ F^jo$\"1+++m2kk>F17$Fcjo$\"1+++*z7>$>F17$Fhjo$\"1+++sDc,>F17$F][p$\"1+ ++9u=t=F17$Fc[p$\"1+++_V[Y=F17$Fh[p$\"1+++Z'=7#=F17$F]\\p$\"1+++SU?(z \"F17$Fb\\p$\"1+++d1HuBas;F17$F`^p$ \"1+++0MAa;F17$Fe^p$\"1+++\"))\\kj\"F17$Fj^p$\"1+++_!z\">;F17$F__p$\"1 +++^JP-;F17$Fd_p$\"1+++h%)*fe\"F17$Fi_p$\"1+++L[-q:F17$F^`p$\"1+++!3DW b\"F17$Fc`p$\"1+++(pu\"R:F17$Fh`p$\"1+++g8DC:F17$F]ap$\"1+++&zM'4:F17$ Fbap$\"1+++)R1`\\\"F17$Fgap$\"1+++q\"\\7[\"F17$F\\bp$\"1+++&QZuY\"F17$ Fabp$\"1+++ul)QX\"F17$Ffbp$\"1+++mKbS9F17$F[cp$\"1+++X]VF9F17$F`cp$\"1 +++J-_99F17$Fecp$\"1+++3!)z,9F17$Fjcp$\"1+++l\"e#*Q\"F17$F_dp$\"1+++<7 *oP\"F17$Fddp$\"1+++H\")ok8F17$Fidp$\"1+++d/k_8F17$F^ep$\"1+++4,uS8F17 $Fcep$\"1+++t&z*G8F17$Fhep$\"1+++9:N<8F17$F]fp$\"1+++7\"\\eI\"F17$Fbfp $\"1+++xcY%H\"F17$Fffp$\"1+++g]>$G\"F17$Fjfp$\"1+++y5.s7F17$F^gp$\"1++ +9\"o4E\"F17$F*$\"1+++&H++D\"F17$Ffgp$\"1+++-/sc7F17$F[hp$\"1+++t4Pj7F 17$F_hp$\"1+++KR&*p7F17$Fdhp$\"1+++j1Zw7F17$Fihp$\"1+++&*G#HG\"F17$F^i p$\"1+++')>J*G\"F17$Fcip$\"1+++C%RcH\"F17$FP$\"1+++2l!>I\"F17$F[jp$\"1 +++aY638F17$F`jp$\"1+++b]E98F17$Fdjp$\"1+++*)*e.K\"F17$Fijp$\"1+++)f(R E8F17$F^[q$\"1+++x?QK8F17$F\\r$\"1+++)[8$Q8F17$Ff[q$\"1+++;H>W8F17$Fd[ l$\"1+++b8-]8F17$F^\\q$\"1+++W)*zb8F17$Fadl$\"1+++.$H:O\"F17$Ff\\q$\"1 +++\"p5sO\"F17$Fh_m$\"1+++%)[%GP\"F17$F^]q$\"1+++)yK%y8F17$Fc]q$\"1+++ E_(RQ\"F17$Fh]q$\"1+++NIZ*Q\"F17$F]^q$\"1+++()p#\\R\"F17$Fb^q$\"1+++') yL+9F17$Fccn$\"1+++ckq09F17$Fj^q$\"1+++WM.69F17$F__q$\"1+++U&>jT\"F17$ Fd_q$\"1+++jac@9F17$Fi_q$\"1+++g=xE9F17$F]`q$\"1+++(RR>V\"F17$Fb`q$\"1 +++!poqV\"F17$Fg`q$\"1+++\"Qg@W\"F17$F\\aq$\"1+++]]@Z9F17$Faaq$\"1+++) HLAX\"F17$Ffaq$\"1+++$o:sX\"F17$F[bq$\"1+++!yi@Y\"F17$F`bq$\"1+++>^2n9 F17$Febq$\"1+++QK&>Z\"F17$Fjbq$\"1+++Swzw9F17$F^cq$\"1+++\\)3;[\"F17$F ccq$\"1+++QtQ'[\"F17$Fhcq$\"1+++,O8\"\\\"F17$F]dq$\"1+++'4[e\\\"F17$Fb dq$\"1+++'HJ0]\"F17$Fgdq$\"1+++RO=0:F17$F\\eq$\"1+++qb!)4:F17$Faeq$\"1 +++1vR9:F17$Ffeq$\"1+++z)f*=:F1Fjeq-F$6$7`qFj[t7$Faho$\"1+++!)QhTFF17$ Fgho$\"1+++D$z(\\GF17$F\\io$\"1+++%3MF17$Fh[p $\"1+++^uEpMF17$F]\\p$\"1+++@8'y^$F17$Fb\\p$\"1+++.bSlNF17$Fg\\p$\"1++ +\\B-7OF17$F\\]p$\"1+++pO\"yl$F17$Fa]p$\"1+++T^'Gq$F17$Ff]p$\"1+++y'\\ su$F17$F[^p$\"1+++Z)H5z$F17$F`^p$\"1+++M)fU$QF17$Fe^p$\"1+++`o)p(QF17$ Fj^p$\"1+++kBD>RF17$F__p$\"1+++UH4hRF17$Fd_p$\"1+++s5a-SF17$Fi_p$\"1++ +-diVSF17$F^`p$\"1+++\"zsV3%F17$Fc`p$\"1+++gc![7%F17$Fh`p$\"1+++&RX\\; %F17$F]ap$\"1+++y5\"[?%F17$Fbap$\"1+++v+UWUF17$Fgap$\"1+++)>)y$G%F17$F \\bp$\"1+++G*HHK%F17$Fabp$\"1+++=&e=O%F17$Ffbp$\"1+++phe+WF17$F[cp$\"1 +++'4C\"RWF17$F`cp$\"1+++wE[xWF17$Fecp$\"1+++%[rc^%F17$Fjcp$\"1+++:%*p `XF17$F_dp$\"1+++4Zd\"f%F17$Fddp$\"1+++X]IHYF17$Fidp$\"1+++kv*om%F17$F ^ep$\"1+++V*eVq%F17$Fcep$\"1+++.apTZF17$Fhep$\"1+++$y7*yZF17$F]fp$\"1+ ++Sl,;[F17$Fbfp$\"1+++/=,`[F17$Fffp$\"1+++$Q.**)[F17$Fjfp$\"1+++5epE\\ F17$F^gp$\"1+++TLRj\\F17$F*Fb^q7$Ffgp$\"1+++o&>l.&F17$F[hp$\"1+++Kc&H2 &F17$F_hp$\"1+++m:J4^F17$Fdhp$\"1+++!e!fX^F17$Fihp$\"1+++%p&z\"=&F17$F ^ip$\"1+++)yHz@&F17$Fcip$\"1+++!e&*RD&F17$FP$\"1+++ic***G&F17$F[jp$\"1 +++%[KfK&F17$F`jp$\"1+++)Q3=O&F17$Fdjp$\"1+++*eDwR&F17$Fijp$\"1*****z? 'QLaF17$F^[q$\"1+++_A4paF17$F\\r$\"1+++\\cu/bF17$Ff[q$\"1+++D#[.a&F17$ Fd[l$\"1+++L?_'F17$F]dq$\"1+++C\"ymb'F17$Fbdq$\"1+++(>JF17$ $!1nm;/R=0vF1$\"1mm;>K'*)>$F17$$!1++]P8#\\4(F1$\"1++]Kd,\"G$F17$$!1nm; /siqmF1$\"1mm;fX(eO$F17$$!1++](y$pZiF1$\"1++]U7Y]MF17$$!1LLL$yaE\"eF1$ \"1LLLV!pu`$F17$$!1nmm\">s%HaF1$\"1mmmhb59OF17$$!1+++]$*4)*\\F1$\"1*** ***H,Q+PF17$$!1+++]_&\\c%F1$\"1+++]*3qy$F17$$!1+++]1aZTF1$\"1+++q=\\qQ F17$$!1nm;/#)[oPF1$\"1mm;fBIYRF17$$!1MLL$=exJ$F1$\"1LLLj$[k.%F17$$!1ML LL2$f$HF1$\"1LLL`Q\"G6%F17$$!1++]PYx\"\\#F1$\"1++]s]k,UF17$$!1MLLL7i)4 #F1$\"1LLL`dF!G%F17$$!1++]P'psm\"F1$\"1++]sgamVF17$$!1++]74_c7F1$\"1++ ]&F17$$\"1)****\\2'HKHF1$\"1*****\\ @fkG&F17$$\"1lmmmZvOLF1$\"1LLL`4Nn`F17$$\"1+++]2goPF1$\"1+++],s`aF17$$ \"1KL$eR<*fTF1$\"1nm;zM)>`&F17$$\"1+++])Hxe%F1$\"1+++qfa(>OF17$F`cv$ \"1nm;>K'*)p$F17$Fecv$\"1++]Kd,\"y$F17$Fjcv$\"1nm;fX(e'QF17$F_dv$\"1++ ]U7Y]RF17$Fddv$\"1MLLV!pu.%F17$Fidv$\"1nmmhb59TF17$F^ev$\"1+++I,Q+UF17 $Fcev$\"1+++]*3qG%F17$Fhev$\"1+++q=\\qVF17$F]fv$\"1nm;fBIYWF17$Fbfv$\" 1MLLj$[k`%F17$Fgfv$\"1LLL`Q\"Gh%F17$F\\gv$\"1++]s]k,ZF17$Fagv$\"1LLL`d F!y%F17$Ffgv$\"1++]sgam[F17$F[hv$\"1++]&F17$F_iv$\"1mm;f`@'G&F17$Fdiv$\"1 ++]nZ)HO&F17$Fiiv$\"1mmmJy*eW&F17$F^jv$\"1+++S^bJbF17$Fcjv$\"1+++0TN:c F17$Fhjv$\"1++]7RV'p&F17$F][w$\"1*****\\@fky&F17$Fb[w$\"1LLL`4NneF17$F g[w$\"1+++],s`fF17$F\\\\w$\"1nm;zM)>.'F17$Fa\\w$\"1+++qfa'F17$F[]w$\"1++]#G2AG'F17$F`]w$\"1LLL$)G[kjF17$Fe]w$\"1++] 7yh]kF17$Fj]w$\"1nmm')fdLlF17$F_^w$\"1nmm,FT=mF17$Fd^w$\"1ML$e#pa-nF17 $Fi^w$\"1,++Sv&)znF17$F^_w$\"1LLLGUYooF17$Fc_w$\"1mmm1^rZpF17$Fh_w$\"1 ++]sI@KqF17$F]`w$\"1++]2%)38rF17$F`eoFfaqFd`w-%+AXESLABELSG6$%!GF[jw-% %VIEWG6$;F(F`eo%(DEFAULTG" 1 2 0 1 0 2 9 1 4 2 1.000000 45.000000 45.000000 0 }}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 142 "To prove rigorou sly that there is a single root between the green lines, and that it i ncreases as a function of c, we compute the derivative " }{XPPEDIT 18 0 "diff(F,lambda);" "6#-%%diffG6$%\"FG%'lambdaG" }{TEXT -1 212 " in an area covering the area between the lines, verify that it is positive, and verify that the characteristic polynomial has different signs on \+ the two lines. All calculations are done in interval arithmetic. " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 1896 "VerifyMaxEvBounds := proc( charpo ly, lowOffset::numeric, highOffset::numeric, slope::numeric, \n \+ csteps::integer, lambdasteps::integer, crange) \n \+ local i,j, deltac, deltalambda, charpolyLowI, charpolyHighI, \n \+ deriv, derivI, curDeriv, curCharHigh, curCharLow, d, failflag, curCra nge, curLambdaRange; \n deltac := (crange[2]-crange[1])/csteps;\n de ltalambda := (highOffset- lowOffset)/lambdasteps;\n charpolyLowI := s ubs( lambda = lowOffset + slope*(c-crange[1]), charpoly);\n charpolyL owI := inapply(charpolyLowI, c );\n charpolyHighI := subs( lambda = h ighOffset + slope*(c-crange[1]), charpoly);\n charpolyHighI := inappl y(charpolyHighI, c );\n deriv := diff( charpoly, lambda ); \n derivI := subs( lambda = d + slope*(c-crange[1]), deriv); \n derivI := inap ply(derivI, [c,d] );\n failflag := false;\n\n for i from 0 to csteps -1 do \n curCrange := map(evalf, [crange[1]+ deltac*i, crange[1]+ d eltac*(i+1)]);\n curCharHigh := charpolyHighI( curCrange);\n cur CharLow := charpolyLowI( curCrange);\n if curCharHigh[1]*curCharH igh[2] <= 0 or\n curCharLow[1]*curCharLow[2] <= 0 or\n c urCharLow[2]*curCharHigh[1] >= 0 then \n failflag := true;\n \+ print(`Sign change test failed for interval`, curCrange,curCharLow,cu rCharHigh); \n fi;\n for j from 0 to lambdasteps-1 do \n cu rLambdaRange := map(evalf, [lowOffset+j*deltalambda+(curCrange[1]-cran ge[1])*slope, \n lowOffset+(j+1) *deltalambda+(curCrange[1]-crange[1])*slope]);\n curDeriv := deri vI( curCrange, curLambdaRange);\n if curDeriv[1]*curDeriv[2] <= 0 or curDeriv[1] < 0 then \n print(`Derivative may be not positi ve for c,lambdain the ranges`, curCrange, curLambdaRange);\n fa ilflag := true; \n fi; \n od;\n if i mod 10 = 0 then print (i);fi; \n od;\n if not failflag then print( `All tests succeeded`); fi;\nend:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 112 "VerifyMaxEvBounds( CharPolyMainFactor, 47/10 0, 52/100, 20/100, \n 160, 10, [0,85/100]); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"!" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#5" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#?" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#I" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#S " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#]" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#g" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#q" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#\"#!)" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#!*" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#\"$+\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"$5\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"$?\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"$I\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"$S\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"$]\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6&%ESign~change~test~failed~for~intervalG7$$\"+++D\"=)!#5$\"++]PM# )F'7$$!+,aM&H#!#7$\"+,+]%3\"!#:7$$\"+*H2<3%F-$\"+,h\\PqF-" }}{PARA 11 "" 1 "" {XPPMATH 20 "6&%ESign~change~test~failed~for~intervalG7$$\"++] PM#)!#5$\"+++](G)F'7$$!+,hq'H#!#7$\"+,+oj8!#97$$\"+**4*H8%F-$\"+,Bf0rF -" }}{PARA 11 "" 1 "" {XPPMATH 20 "6&%ESign~change~test~failed~for~int ervalG7$$\"+++](G)!#5$\"++]iS$)F'7$$!+,-\"zH#!#7$\"+,+'3k#!#97$$\"+*fb Z=%F-$\"+,2CurF-" }}{PARA 11 "" 1 "" {XPPMATH 20 "6&%ESign~change~test ~failed~for~intervalG7$$\"++]iS$)!#5$\"+++v$R)F'7$$!+,h&*)H#!#7$\"+,+A SR!#97$$\"+*f.qB%F-$\"+,XWVsF-" }}{PARA 11 "" 1 "" {XPPMATH 20 "6&%ESi gn~change~test~failed~for~intervalG7$$\"+++v$R)!#5$\"++](oW)F'7$$!+,;% )*H#!#7$\"+,+'>E&!#97$$\"+*zQ(*G%F-$\"+,u?8tF-" }}{PARA 11 "" 1 "" {XPPMATH 20 "6&%ESign~change~test~failed~for~intervalG7$$\"++](oW)!#5$ \"+++++&)F'7$$!+,\\c+B!#7$\"+,+G1m!#97$$\"+*>kHM%F-$\"+,F`$Q(F-" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "VerifyMaxEvBounds( CharPoly MainFactor, 47/100+17/100, 52/100+17/100, 20/100, \n \+ 50, 10, [85/100,1]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" } }{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"!" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#5" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#?" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#\"#I" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#S" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#%4All~tests~succeededG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" } }}{SECT 0 {PARA 0 "" 0 "" {TEXT -1 83 "(3). We \"deflate\" the chara cteristic polynomial (that is, divide by the monomial " }{XPPEDIT 18 0 "lambda - mu" "6#,&%'lambdaG\"\"\"%#muG!\"\"" }{TEXT -1 30 ") in th e symbolic form, with " }{XPPEDIT 18 0 "mu" "6#%#muG" }{TEXT -1 5 " an d " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 53 " as the indeterminates. Next, we verify that for all " }{XPPEDIT 18 0 "c = 0..1" "6#/%\"cG;\" \"!\"\"\"" }{TEXT -1 14 ", and for all " }{XPPEDIT 18 0 "lambda = 0.5. .0.613" "6#/%'lambdaG;$\"\"&!\"\"$\"$8'!\"$" }{TEXT -1 146 ", all root s of the deflated polynomial are inside the circle of radius 0.5 cente red at 0 in the complex plane, that is, have magnitudes less than " } {XPPEDIT 18 0 "mu(c)" "6#-%#muG6#%\"cG" }{TEXT -1 10 " for any " } {XPPEDIT 18 0 "c > 0" "6#2\"\"!%\"cG" }{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "Symbolic deflat ion; if we substitute a pair " }{XPPEDIT 18 0 "c, mu(c)" "6$%\"cG-%#mu G6#F#" }{TEXT -1 89 " we get the deflated polynomial for a specific va lue of c. Normalize deflated polynomial." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 92 "deflatedCharPoly := collect( expand( deflate(CharPoly MainFactor, lambda, mu)), lambda)/64^5;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%1deflatedCharPolyG,B*$)%'lambdaG\"\"%\"\"\"#\"\"\"\"+C=ut5*&,( #!#D\"#KF,%\"cG#!\"&F2%#muGF,F,)F(\"\"$F*F+*&,.F3#\"#<\"$7&#\"#@\"$G\" F,*$)F3\"\"#F*#!\"\"F=F6F0*&F6F,F3F,F4*$)F6FCF*F,F,)F(FCF*F+*&,4FA#\" \"(\"%#>)#F1\"%[?F,F3#!#8FNFFF;F6F>*&F6F*FBF*FDFGF0*&FHF*F3F*F4*$)F6F8 F*F,F,F(F,F+F3#FE\"/K))3sV=N*&FHF*FBF*#FE\"-))Q\"ev\\&#F,\"./6^Y!)R%F, F6#F1\"._bDB!*>#*$)F6F)F*F+FS#FM\".3A-$4'z)*&FVF*F3F*#F5\",o$Q(fV$FAFW FF#FRF]oFT#F " 0 "" {MPLTEXT 1 0 1248 "TestDeflated := proc(csteps::integer, lambdasteps::integer, hig hOffset, lowOffset, slope, crange ) \n local cf,i,j,m, MJ, cinterv, d eflatedinterv, lambdainterv, deltac, deltalambda; \n deltac := (cran ge[2]-crange[1])/csteps;\n deltalambda := (highOffset- lowOffset)/lam bdasteps;\n for i from 0 to 3 do \n cf[i] := inapply( coeff(defla tedCharPoly,lambda,i), c, mu);\n od;\n for i from 0 to csteps-1 do\n cinterv := [ crange[1]+i*deltac, crange[1]+(i+1)*deltac];\n for j from 0 to lambdasteps-1 do \n lambdainterv := [ lowOffset + (c interv[1]-crange[1])*slope + j*deltalambda,\n l owOffset + (cinterv[2]-crange[1])*slope + (j+1)*deltalambda];\n d eflatedinterv := 0;\n for m from 0 to 2 do \n deflatedinte rv := deflatedinterv + eval(cf[m](cinterv,lambdainterv))*lambda^m; \+ \n od;\n # fix for cf[3]: stupid inapply does not always ret urn an interval\n deflatedinterv := deflatedinterv + [1.0,1.0]*la mbda^3;\n MJ := IntervMardenJury( deflatedinterv, lambda, 0.5);\n if op(2, MJ[1]) > 0 or op(1, MJ[2]) < 0 or op(1,MJ[3]) < 0 \n \+ then ERROR(`test failed for interval `, cinterv); \n fi;\n \+ od;\n if i mod 10 = 0 then print(i); fi;\n od;\n print(`All t ests passed`);\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 22 "This may \+ take a while " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "TestDeflated(10,10 ,0.5,0.47, 0.2, [0.0,1.0] );" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"! " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%1All~tests~passedG" }}}}{SECT 0 {PARA 0 "" 0 "" {TEXT -1 19 "Special case: k = 3" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 99 "All we have to do is to run Marden-Jury test for the deflated characteristic polynomial for k = 3" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "mu3 := fsolve( subs( c = cos(2*Pi/3), CharPolyMainFac tor ), lambda, 0.3..1); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$mu3G$\" +1diVS!#5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 19 "Verify precision: \+ " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "inapply( subs( c = cos(2*Pi/3), CharPolyMainFactor ), lambda)(mu3 + 1e-7);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$$\"+******H5!#=$\"+,++q6F&" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "inapply( subs( c = cos(2*Pi/3),CharPolyMainFactor \+ ), lambda)(mu3 - 1e-7);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$$!+,++!= \"!#=$!+******H5F&" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "inter vmu3 := [mu3 - 1e-7, mu3 + 1e-7];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#> %*intervmu3G7$$\"+1ZiVS!#5$\"+1niVSF(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "Construct the interval deflated polynomial" }{MPLTEXT 1 0 313 " \n for i from 0 \+ to 2 do \n cf3[i] := inapply( coeff(subs( c = cos(2*Pi/3), deflate (CharPolyMainFactor/64^5 ,lambda,mu)),lambda,i), mu);\n od: \n cf3[ 3] := mu -> [1.0,1.0]: \n deflatedinterv3 := 0:\n for i from 0 to 3 \+ do \n deflatedinterv3:= deflatedinterv3 + eval(cf3[i](intervmu3)) *lambda^i; \n od:\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 93 "Verify that all other roots have magnitude < 0.3: the Marden Jury test for r adius 0.3 passes:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "IntervMardenJu ry(deflatedinterv3, lambda, 0.3);" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#- %&TABLEG6#7%/\"\"\"7$$!+0+++5!\"*$!+x********!#5/\"\"#7$$\"+_********F /$\"+7+++5F,/\"\"$7$$\"+-********F/$\"+E+++5F," }}}}}{SECT 0 {PARA 4 " " 0 "" {TEXT -1 67 "Calculation of the largest eigenvalues with guar anteed precision " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 370 "This functio n produces a table of approximate values of the eigenvalue with given \+ precision for use with interval arithmetics in the C part of the analy sis code; to avoid conversion problems, we \nwrite intervals as 3 long integers: low, high, scale, which represent the interval [low/scale, \+ high/scale].\n The first value is the limit value for infinity (comp uted with " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 319 " set to 1 in t he char. polynomial). The result is a C function written to a file; i f the file name is `default`, then the output is written to the standa rd output.\nThe argument of the function is valence, the function retu rns the interval value for the largest eigenvalue. The body just cont ains a big static array. \n" }{MPLTEXT 1 0 2243 "\nComputeEigenvalues \+ := proc(N::integer, eps::numeric, fname) \n local K, intervCharPoly, \+ intervPi, intervc, \n expandedCharPoly, approxEV, r, deflatedCharPol y, i,marTable, cK; \n global SimpleChaprPolyScaled; \n Digits := 15 ;\n intervCharPoly := inapply( CharPolyMainFactor ,lambda, c);\n # u se arccos to make sure we get an interval for Pi\n intervPi := Interv al_times([2.0,2.0],Interval_arccos([0,0])); \n intervc := K -> Inte rval_cos(Interval_times( [2*intervPi[1], 2*intervPi[2]], Interval_reci procal(map(evalf,[K,K])))); \n fprintf(fname, `virtual Float Eigenva lue(int K) \{\\n`);\n fprintf(fname, ` static INTEGER64 EV[] = \{\\n` );\n # do infinty\n expandedCharPoly := subs( c = 1, CharPolyMainFa ctor);\n approxEV := fsolve( expandedCharPoly, lambda, lambda=0.5..1) ;\n if op(2, intervCharPoly(approxEV-eps, 1)) > 0 \n or op(1,inter vCharPoly(approxEV+eps, 1)) < 0 then \n ERROR(`fsolve precision fai lure for infinity`);\n fi;\n fprintf( fname, `CONST64(%d),CONST64(%d ),CONST64(%d),\\n`, op(1,approxEV-eps), op(1,approxEV+eps),10^(-op(2,a pproxEV)));\n # skip 1,2\n fprintf( fname, `CONST64(0),CONST64(0),CO NST64(0), CONST64(0),CONST64(0),CONST64(0),\\n`);\n for K from 3 to N do \n if (K - 3) mod 100 = 0 then print(K);\n fi; \n cK := intervc(K);\n expandedCharPoly := subs( c = cos(2*Pi/K), CharPo lyMainFactor);\n approxEV := fsolve( expandedCharPoly, lambda, lamb da=0.25..1); \n # check that the precision is at least eps \n if op(2, intervCharPoly(approxEV-eps, cK)) > 0 \n or op(1, interv CharPoly(approxEV+eps, cK)) < 0 then \n ERROR(`fsolve precision \+ failure for K =`, K);\n fi;\n # writing out f.p. numbers means r elying on the parser of \n # the compiler, which will probably do r ounding incorrectly when converting to the \n # binary format; inst ead, we write out integers and assume that the interval numbers are\n \+ # intialized correctly from integers, i.e. conversion is done by th e processor in the correct \n # rounding mode;\n fprintf(fname, \+ `CONST64(%d),CONST64(%d),CONST64(%d),\\n`, op(1,eval(approxEV-eps)), o p(1,eval(approxEV+eps)),10^(-op(2,approxEV)));\n od;\n fprintf( fname, `CONST64(0)\};\\n return Float(EV[3*K],EV[3*K+1])/Float(EV[3*K+ 2]) ;\\n\}\\n\\n`);\n NULL;\nend: " }}}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 72 "The derivative of the largest eigenvalue with respect to \+ c at infinity. " }}{PARA 0 "" 0 "" {TEXT -1 190 "To establish C1-conti nuity for all valences, we need to analyze behavior of the magnitude o f the largest eigenvalue as the function of the valence, as the val ence increases to infinity ( " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 53 " approaches 1). We estimate a constant R, such that " }{XPPEDIT 18 0 "abs(lambda-lambda[infinity]) < R*abs(c-1);" "6#2-%$absG6#,&%'lam bdaG\"\"\"&F(6#%)infinityG!\"\"*&%\"RGF)-F%6#,&%\"cGF)\"\"\"F-F)" } {TEXT -1 77 ", sufficiently close to 1. This constant can be taken to be the maximum of " }{XPPEDIT 18 0 "abs(diff(lambda,c))" "6#-%$absG6 #-%%diffG6$%'lambdaG%\"cG" }{TEXT -1 36 ", or, equivalently, as maxim um of " }{XPPEDIT 18 0 "abs( diff(c,lambda))^(-1)" "6#)-%$absG6#-%%di ffG6$%\"cG%'lambdaG,$\"\"\"!\"\"" }{TEXT -1 52 "; as the characterist ic polynomial is quadratic in " }{XPPEDIT 18 0 "c" "6#%\"cG" }{TEXT -1 50 ", the latter is relatively easy to compute. Once " }{XPPEDIT 18 0 "B" "6#%\"BG" }{TEXT -1 36 " is known, we can estimate the size \+ " }{XPPEDIT 18 0 "epsilon" "6#%(epsilonG" }{TEXT -1 22 " of the inter val for " }{XPPEDIT 18 0 "lambda" "6#%'lambdaG" }{TEXT -1 7 " near " }{XPPEDIT 18 0 "lambda[infinity]" "6#&%'lambdaG6#%)infinityG" }{TEXT -1 134 " , such that if the characteristic map is injective and regula r for all these values, it is sufficient to establish C1-continuity fo r " }{XPPEDIT 18 0 "K > K[0]" "6#2&%\"KG6#\"\"!F%" }{TEXT -1 8 ", wher e " }{XPPEDIT 18 0 "epsilon/R < 1-cos(2*Pi/K[0]);" "6#2*&%(epsilonG\" \"\"%\"RG!\"\",&\"\"\"F&-%$cosG6#*(\"\"#F&%#PiGF&&%\"KG6#\"\"!F(F(" } {TEXT -1 1 "." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "lambdac := \+ [];\nintervlambdac := inapply( -Fc/Flambda, [c,lambda]):" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%(lambdacG7\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 70 "This shows that the maximum of the derivative can be estimated \+ by 0.5:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 253 " for i from 0 to 19 do \n for j from 0 to 9 do \n cinterv := [0.9+i*0.005, 0.905+i*0 .005]; \n lambdac := Interval_union(lambdac,intervlambdac( cinter v, [0.47+0.2*cinterv[1]+j*0.005, 0.48+0.2*cinterv[2]+j*0.005]));\n \+ od;\n od; eval(lambdac);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$$\"+wpa /&*!#6$\"+2*4M1%!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "R := op(2, lambdac);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"RG$\"+2*4M1%!# 5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 15 "Code generation" }}{PARA 0 "" 0 "" {TEXT -1 53 "A class with three memeber functions are generated:\n " }{TEXT 256 23 " Float Eigenvalue(int K)" }{TEXT -1 28 " computes the eigenvalues, \n" }{TEXT 257 58 "void EigenvectorReal(Float c, Float lambda, Float EvRe[ ]) " }{TEXT -1 68 "initializes an array for the real part of the compl ex eigenvector, " }{TEXT 258 62 "void EigenvectorImaginary(Float c, F loat lambda, Float EvIm[])" }{TEXT -1 108 " initializes the array for \+ the complex part. \nMemory for arrays should be allocated by the calli ng function." }}{PARA 0 "" 0 "" {TEXT -1 267 "The output is written to a file; if the name is `default`, it is written to the standard outp ut (warning: for some reason, writing to standard output is terribly s low; writing to a file and then looking at it in an editor is much mo re efficient. All functions use " }{TEXT 259 5 "Float" }{TEXT -1 158 " as the name of the class for the interval numbers.\nIt is assumed t o have explicit casts from 64-bit integers, standard arithmetics ope rations, and macros " }{TEXT 260 3 "FR " }{TEXT -1 6 " and " }{TEXT 261 4 "Fdiv" }{TEXT -1 9 ", (see " }{MPLTEXT 0 21 14 "ConvertToFloat " }{TEXT -1 15 " for details).\n" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "OutputFile := \"velho.cpp\":" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 280 "Generate the header and the constructor for the class; specify the name of the class, the inner and outer diameters of the layer for the char. map, the radius of the neighborhood on wehich the eigenvect or is defined, and the name of the corresponding scheme\non the regula r complex" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }{MPLTEXT 1 0 57 "MakeClas sHeader( OutputFile, \"Velho\", 2,4,3, \"EightBox\"):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "ComputeEigenvalues(450, 1e-10, Outp utFile):" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"$" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"$.\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"$.#" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#\"$.$" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"$.%" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 72 "This function is need ed for computing the eigenvectors \"near infinity\"." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 269 "fprintf( OutputFile, `virtual void EigenvalueRa nge(Float c, Float& lambdamin, Float& lambdamax) \{\\n`):\nfprintf( Ou tputFile, `lambdamax = Eigenvalue(0); lambdamin = Eigenvalue(0) -\n ex actfloor(FR(CONST64(%d), CONST64(%d))*(Float(1)-c));\\n \}\\n`, op(1, R), 10^(-op(2,R))):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "eval m( subs( \{ c= 0, s = 1, lambda = 1/2\}, eval(EV) ));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'vectorG6#7.\"\"\",&F'F'%\"IGF'\"\"#,&F*F'F)F',& F*F'F)F*,&F'F'F)F*\"\"$,&F.F'F)F',&F.F'F)F*,&F.F'F)F.,&F*F'F)F.,&F'F'F )F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "GenerateEigenvectorC ode(EV, OutputFile);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "End of th e class " }{MPLTEXT 1 0 26 "fprintf(OutputFile, `\};`):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "fclose(OutputFile):" }}}}}{MARK "7 \+ 5 3 5 2 0" 0 }{VIEWOPTS 1 1 0 3 4 1802 }