|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| Full Text | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
PAGE 1 SIMULATINGELECTROMAGNETICWAVEPROPAGATIONUSINGMULTI-SCALE TIME-DOMAINMETHODSANDITERATIVEALGORITHMS BY KERRYKEYS AThesis SubmittedtotheDivisionofNaturalSciences NewCollegeofFlorida inpartialfulllmentoftherequirementsforthedegree BachelorofArts UnderthesponsorshipofGeorgeRuppeiner Sarasota,Florida May,2011 PAGE 2 Contents 1Introduction 1 1.1GoalsofNumericalSimulation.............................1 1.2HistoryandApplicationofFinite-DierenceTime-DomainMethods.........2 2Wavelets 4 2.1ABriefOverviewofFourierAnalysis..........................4 2.2RepresentingFunctionswithaMultiresolutionBasis..................7 2.3WaveletsasaFiner-scaleCompaniontoScalingFunctions..............9 2.4FamiliesofWavelets...................................12 2.5ScalingCoecients...................................15 3RepresentingFunctionsandDerivativeOperatorswithaWaveletBasis20 3.1ProjectingFunctionsintoaWaveletBasis.......................20 3.2ContinuousandDiscreteWaveletMoments......................23 3.3ApproximatingProjectionCoecientsusingNumericalQuadrature..........25 3.4Finite-DierenceDerivativeOperators.........................26 3.5ProjectingDerivativeOperatorsintoWaveletSpace..................29 4Maxwell'sEquationsinaSchrodinger-likeForm31 4.1Maxwell'sEquations...................................31 4.2ConstructingaHamiltonianOperator..........................32 i PAGE 3 5IterativeAlgorithms34 5.1PowerIteration.....................................35 5.2KrylovSubspaceTechniques..............................37 5.3TheLanczosAlgorithm.................................38 6ExampleSimulation40 7Conclusion 45 ASizeRequirementonScalingCoecients46 ii PAGE 4 Abstract SIMULATINGELECTROMAGNETICWAVEPROPAGATIONUSINGMULTI-SCALE TIME-DOMAINMETHODSANDITERATIVEALGORITHMS KerryKeys NewCollegeofFlorida,2011 Techniqueswhichimprovetheeciencyofnumericalsimulationsusedtoevaluatepartialdierential equationsareespeciallyimportantincreatingrealisticmodelsofphysicalphenomena.Inthisthesis Iexaminetheimplementationofawaveletbasisinsimulatingthepropagationofelectromagnetic wavesthroughadielectricmedium.Ireviewthepropertiesofawaveletbasisanditsadvantages overthetraditionalFourierbasis.Ithenexplainthegeneralmethodofprojectingfunctionsand operatorsintoagivenwaveletbasis,andfollowtheconstructionofaHamiltonianoperatorwhich denesthebehaviorofelectromagneticeldsinagivensystem.Thismaterialsupportstheexample propagationofanelectromagneticwavethroughadielectricmaterialusingavariationofaniterative computationalmethodknownastheLanczosAlgorithm. SponsoredbyGeorgeRuppeiner NewCollegeofFloridaDivisionofNaturalSciences PAGE 5 Chapter1 Introduction 1.1GoalsofNumericalSimulation AnalyticsolutionstoMaxwell'sequationsaregenerallylimitedtoproblemsinvolvingspecialcases ofsymmetry,sosolving'real-world'problemsincludingthosethatinvolveasymmetricorirregular materialgeometriesnecessitatesnumericalsimulation.Developingnewnumericaltechniquesallows ustoworkonavarietyofproblemsmoreecientlyandtoincreasetheclassofproblemswhich aretractableusingavailablehardwareresources.However,betternumericalmodelsaren'teasyto comeby-currentchallengesincludedevelopingmodelsthatcanpropagateelectromagneticwaves throughlossymedia,acrossdiscontinuousdielectricfunctionsaswhenwemovefromonemedium toanother,andinthepresenceofcurrentandchargesources. WecanbegintoaddresstheseissuesbyextendingthepopularFinite-DierenceTime-Domain FDTDmethodtousemulti-scalebasisfunctionsthatis,functionswhichcanberepresentedatone time-orspace-scaleasalinearcombinationofthesamebasisfunctionsatanerscaleandsimulating thepropagationofelectromagneticwavesusingavariationofaniterativecomputationalmethod knownastheLanczosAlgorithm.Byimprovinguponthesemethodswecanhopetousereadilyavailablecomputationalresourcestomodeldemandingproblemslikethescatteringandabsorption 1 PAGE 6 oflightbymetalnanoparticles.Doingsomayenableresearcherstobetterunderstandandpredict expectedvaluesforsurfaceplasmoneectsandRamanSpectroscopy,anditmaybeaboonto nanophotonicsresearchingeneral[1]. 1.2HistoryandApplicationofFinite-DierenceTime-Domain Methods TheFinite-DierenceTime-DomainFDTDmethodwaslargelymotivatedbyanddevelopedon behalfofcold-wareramilitaryconcerns.Importantapplicationsincludeddevelopingstealthbombers, improvingmissileradar-guidancesystems,andmodellinghowtheelectromagneticpulseEMPfrom anuclearbombaectselectroniccomponents[2].Computationalmodellingwasadvantageousin theseprojectsas,forexample,simulatinghowanEMPaectselectricalcomponentsallowedmilitary plannerstooutlinehowanuclearattackmightaectelectricalsystemsinanattackcarriedout eitherbyoragainsttheUnitedStates,whileengineerswereabletodesignEMP-tolerantelectronics systemsforcriticalapplications.Inthecaseoftheradar-guidancesystems,FDTDhelpedmodelthe slightinuenceofthemissilenoseontheelectromagneticwavesenroutetothemissile'sinternally housedguidancesystems.Bymodellingtheseinteractionsbeforemissilefabrication,manufacturers couldprogramtheguidancesystemtotakeintoaccountthedeectionofthewavesbythenosein ordertomoreaccuratelyndtheirtargetwhilein-ight. FDTDhasbecomeincreasinglyimportantincircuitdesign,asittakesintoaccountonlythe electricandmagneticeldswhensimulatinghowacircuitfunctions.Asitdoesnotneedtoconsider relationsbetweenvoltage,current,andresistance,theFDTDmethodcanbeespeciallyusefulwhen modellingnonlinearactivedevicessuchasdiodes,transistors,orlogicgates[2].FDTDhasalso provenindispensableasthedesignofcircuitsinvolvesever-increasingoperationspeedsanddecreasing chipsize,leadingtocomplicationsinvolvingelectromagneticwavecouplingbetweensignalpaths, radiation,andground-currentartifacts[2].Inadditiontohavingarangeofapplicationsincircuit 2 PAGE 7 design,FDTDalsoperformsaccurately:computedcircuitquantities"typicallyagreewithtextbook valuesontheorderof1%orbetter". AlsoofmajorimportanceistheapplicationofFDTDmethodstoproblemsinnonlinearoptics. Otherattemptstomodelnonlinearopticssometimesincludea"nonphysicalscalarcarrier-envelope function":here,FDTDshowsitsstrengthbysimplypluggingthroughalargenumberofcomputations for E and B ,doingawaywiththenonphysicalsimplication[2].ThuswecanuseFDTDtomodel thevariationofamaterial'sindexofrefractionwithfrequencyatlowopticalpowerlevelslinear dispersionaswellasthevariationofthefrequency-dependentrefractiveindexwithopticalpower nonlineardispersion.TheapplicationofFDTDtomodellingnonlinearcomponentslendsitselfto designingall-opticalswitchingcircuits,andthusitmayeventuallyplayanimportantroleindesigning opticalcomputers. FDTDhasalsofounduseindesigningphotonicband-gapmicrocavitylasersandinapplications wherethepropagationandabsorptionofE/Mwavesthroughamediumbecomeoverlycomplicated. Foranexampleofthelatter,FDTDhasbeenusedtomodelhowcell-phoneradiationisabsorbed byboneandtissueinpeople'sheads,[2]andscaleswelltoothercongurationsofinhomogeneous media. 3 PAGE 8 Chapter2 Wavelets 2.1ABriefOverviewofFourierAnalysis Ofthealready-establishedmethodsusedtomodeltime-varyingelectromagneticeldproblems,the Finite-DierenceTime-DomainFDTDmethodhasseenthemostwidespreaduseinbothscientic andengineeringcircles.Ingeneral,waveletsplayananalogousroleinMulti-ScaleTime-Domain methodsastheroleFourieranalysisplaysinFDTDmethods,soitishelpfultounderstandthe simplercaserstFourieranalysistohaveanideaofwherewe'reheaded.Thebasicgoalof Fourieranalysissimilarlyforwaveletsistorepresentphysicalphenomenaasaweightedsumof basisfunctions.Thatis,givensomefunctionoftime[orspace], f x ,weexpresstheunderlying informationasasumoffunctionsinthefrequency[orwavenumber]domain. InFourierAnalysis,werepresentanarbitraryfunction, f x ,asasumofsineandcosinebasis functions: f x = a 0 2 + 1 X n =1 [ a n cos nx + b n sin nx ] .1 4 PAGE 9 wherethecoecients a n b n aredeterminedby: a n = 1 Z 2 0 f x cos nx d x; .2a b n = 1 Z 2 0 f x sin nx d x: .2b Ifwecanndasuitablealternatebasisandcoecients,thenwecangeneralizetheaboveequations intoaform f = P n N c n f n Onereasontosearchforanalternativebasisisthepossibilityofndingbasisfunctionswhichare simultaneouslylocalizedintime[orspace]andfrequency[wavenumber].TheFourierbasisfunctions arenotlocalizedatallinthetimedomainthatis,sineandcosinewavesaredenedforalltimes f PAGE 10 .[3] Assuch,manysuperimposedsineorcosinewavesarerequiredtoaccuratelyrepresentfunctions whichdonotmaintaintheirbehavioroverlongperiodssay,forexample,apianonoteplayedfora veryshortduration.Forsuchapplicationsabasiswhichislocalizedinboththetime[space]and frequency[wavenumber]domainsprovidesasparserrepresentationofthedesiredfunctionsandthus requirelessstorageandprocessingforsimilarresults. Manywaveletsaredenedtohavecompactsupport{thatis,theyareidenticallyzerooutsidea closed,boundeddomaineitherintimeorspace,thelengthofwhichisknownasthesupportlength ofthewavelet.Assuch,theyaresomewhatlocalizedintime[space]andthusarewell-suitedfor modellingthebehavioroffunctionswhichvarysharplyoverashortspan.Suchjumpsoccurfor exampleinthesimulatedtransmissionofanelectromagneticwaveasitcrossesaboundaryfromone dielectricmaterialtoanother. 6 PAGE 11 2.2RepresentingFunctionswithaMultiresolutionBasis Inpart,FDTDhasbeensosuccessfulbecauseitssimplicityallowsnumericalmodelerstoscaleits usetoavarietyofgeometriesbysimplyreningthespatialscaleatwhichsolutionstotheeld equationsarecalculated.Initstraditionalframework,however,FDTDusesauniformspatialgrid. Thismeansthatateveryspatialpointinthesimulationthegridspacingmustbeassmallasis requiredforaccuracyatthepartofthemodelwhichismostsensitivetochangesinthegridsize. ThusitiscommoninFDTDthatonesensitivedomainofanumericalmodelwillincreasethespatial gridresolutionoverlargeswathsofspacewhichwouldotherwisesubmittoresolutionatacoarser scale. Increasingtheresolutionoftheproblemgridlikewiseincreasesthecomputationalcomplexityof asimulation,sowewouldbebetteroifwecoulduseamethodwhichemploysacoarserspatial gridinareaswherethedielectricfunctiondoesn'tvarysignicantly,usinganerscaleonlywhereit isrequired.Thisisthebasicpremiseofmultiscaleormultiresolutiontimedomaintechniques:we wanttoexpressafunctionatcoarserscalesasalinearcombinationofthesamefunctionprojected intoanerscalespace.Weaccomplishthisbyswitchingourbasisfunctionsfromthefamiliar sine and cosine tofunctionswhichobeyrelationsbetweencoarserandnerscales: j x = X k c k j +1 x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k : M1 M1tellsusthatbycombiningthener-scalebasisfunctionsastheyaretranslatedalongthexaxisinincrementsof k ,wecanrecoverthesamebasisfunctionatthelargerscale.Notethat incrementingjby1increasesthegridresolutionbydoublingthenumberofwaveletsrequiredtospan thesameamountofspace. Thefunction x iscalledthescalingfunction,anditconstitutesthebasicbuildingblockof amultiresolutionbasis.Section2.4describesseveralscalingfunctions,althoughthesimplestisthe scalingfunctionfortheHaarwaveletbasis: 7 PAGE 12 [4]. Thescalingcoecients, c k ,describehowtotranslatethescalingfunctionbetweencoarserandner scales.Infact,thescalingfunctionisitselfuniquelydeterminedbytheseweights. Sincescalingfunctionsconstituteacompletebasisfor L 2 ,[5]wecanrepresentanyarbitrary square-integrablefunctionintermsof x : f x = X j;k j;k j;k x .3 f Notationalconvention g : j;k x =2 j 2 j x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k wheretheprojectioncoecients, j;k ,give j;k x theproperweightatposition 2 j x )]TJ/F11 9.9626 Tf 7.924 0 Td [(k torepresent fx. 8 PAGE 13 2.3WaveletsasaFiner-scaleCompaniontoScalingFunctions Whilewecanrepresentanysquare-integrablefunctionusinganinnitecombinationofweighted scalingfunctions,therealstrengthofwaveletbasesresultsfrombeingabletodescribethe dierence betweenthescalingfunctionsatcoarserandnerscalesintermsofaseparatefunction, x in somecasesthisfunctionisreferredtoasthe"motherwavelet"andthescalingfunctioniscalledthe "fatherwavelet"{hereIwilladopttheconventionofsimplynaming"wavelets"thefunctionswhich spanthedierencebetweenadjacentscalesofthescalingfunction. Thewaveletbasisthenallowsustorelatevaluesgeneratedbyanersamplingofaparticular areainaproblemdomaintothevaluesgeneratedgivenacoarsersamplingofthewholedomain. Thiscanbeespeciallyhelpfulinproblemswhereinformationismoredenseincertain,well-dened areas,asintherepresentationofachirpfunction: 9 PAGE 14 TheHaarscalingfunctionandwaveletsatdierentscales: 10 PAGE 15 [6] Waveletsalsoobeyamulti-scalerelationsimilartoM1: j x = X k d k j +1 x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k : M2 wherethe d k coecientsdeneaconstructionwhichisorthogonaltothescalingfunction,whichat thesametimespansthedierenceindetailcoveredbythescalingfunctionbetweenadjacentscales. Combiningwaveletsandscalingfunctionsallowsustorepresentagivenfunctionasasumofthe two: f x = X k j =0 ;k j =0 ;k x + X j;k j;k j;k x .4 Equations.3and.4bothclaimtorepresentthefunctionfx:eithergivesanexactrepresentationinthelimit j !1 Thedevelopmentofarepresentationintermsofamultiresolutionbasisasoutlinedaboveis particularlyimportant,soI'llemphasizethesepoints:Thedeningpropertyofamultiresolutionbasis isthecharacterizationofthebasisfunctionsbyformulaewhichallowustowritetheexpansionof aphysicalfunctionasalinearcombinationofthesamebasisfunctionsatadierentscale.Using 11 PAGE 16 aseparatefunctiontorepresentthedierencesbetweenonescaleandanothermeansthatwecan computevaluesforamodelatanescaleandthenonlystorethevalueswhichgivethedierence betweenthenestscaleandtheinformationatacoarserscale.Thispropertyiswhatallowsusto partitionthedomainofasimulationandreservethener-scaleapproximationsforthemostsensitive areas.This,inturn,allowsforsimulationswhichavoidthecomputationalcostsassociatedwith calculatingne-scaleapproximationsovertheirentiredomain. 2.4FamiliesofWavelets Thecentralcharacteristicofanymultiresolutionbasisisthatwecanwritearepresentationofa functionatagivenscaleasalinearcombinationofitsrepresentationatnerscales.However,since thereareinnitelymanybasisfunctionswhichsatisfythissimplerequirement,choosingawell-suited basisforaparticularprobleminvolvesunderstandinghowotherpropertiesofpossiblebasisfunctions inuenceissueslikethepresenceofvisualartifactsinwavelet-basedimageprocessing[7],orthe rateofconvergenceofacommonalgorithmusedtoconstructthescalingfunctionfromthescaling coecients[8]. Historically,AlfredHaar'seortscirca1909todescribethespaceofsquare-integrablefunctionsledhimtodevelopasimple'step-function'basiswhichwasthersttoobeythemulti-scale relationswhichcharacterizeallwaveletsandscalingfunctions[9].WhiletheHaarwaveletissimple initsconstruction,isorthonormalandhascompactsupport,thewaveletandscalingfunctionsare discontinuous,whichmakesthemapoorchoiceforrepresentingcontinuousfunctions. ThenextstepinwavelettheorycamewithIngridDaubechies'sdevelopmentofcompactsupport orthonormalwaveletswithvariableregularityinthelate1980's.TheDaubechiesWaveletsareafamily oforthonormalwaveletswhoseregularityandsupportlengthcorrespondtothenumberofconditions whichdenethescalingcoecientsusedtoconstructthescalingfunctionandwavelet.Thesupport lengthofthewaveletandscalingfunctionincreaseswiththenumberofscalingcoecients,sothat thesupportlengthofDaubechieswaveletsisonelessthanthenumberofconditionsimposedonthe 12 PAGE 17 basisthatis,aDaubechieswaveletwithLscalingcoecientswillhavesupportlength=L-1.The numberof vanishingmoments forthewaveletalsoincreaseswiththenumberofconditionsLthat denethescalingcoecients.Thepresenceofvanishingmomentsforthewaveletimpliesanincrease intheabilitytoexactlyrepresentpolynomialfunctionsintermsofthescalingfunctionbasis,andso higher-orderDaubechieswaveletscanbeusedtomoreaccuratelyrepresentirregularfunctionsatthe costofasmallincreaseinthecomputationalresourcesrequiredtocomputethescalingfunctionand wavelet,aswellasaslightincreasesinstoragerequirementsforapplications. [10] [11] 13 PAGE 18 [12] Otherwaveletsdevelopedsincethe1980'shavefoundspecialnichesbyreducingtheasymmetry ofthebasiswithapplicationsinimageprocessingorbyimposingconditionsforvanishingmoments inthescalingfunctionbasisaswellallowingprojectioncoecientstobeapproximatedbysampling fxforcertainvaluesofx,[7]greatlysimplifyingtheimplementationofthesebasisfunctionsin applications. SymmletL=6ScalingFunctionandWavelet.Thesymmletwaveletsandscalingfunctionsareless asymmetricthanthosebelongingtootherwaveletfamilies.[13] Stillotherwaveletsmaybechosenwhichhaveexplicitfunctiondescriptions,simplifyingtheirexpressionandimplementation. 14 PAGE 19 Theabilitytoconstructbiorthogonalwaveletbasescanfurthercomplicateachoiceofbasis.That is,afunctionbasisforwhichwehavetwosetsofwaveletsandtwosetsofscalingfunctions,each satisfying: f x = X k ~ j k j x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k + X k ~ j k j x )]TJ/F11 9.9626 Tf 9.962 0 Td [(k ; .5a = X k j k ~ j x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k + X k j k ~ j x )]TJ/F11 9.9626 Tf 9.962 0 Td [(k ; .5b wheretheprojectioncoecientsaregivenby: j;k = PAGE 20 ourbasisandthentranslatethesepropertiesintoconstraintsonthescalingcoecientsgivenbythe multi-scalerelationsM1andM2. Forexample,wecanconsidertheDaubechieswaveletwithL=4scalingcoecients.Ifwetake j=0tobeourstartingcoarsestresolutionthentheformulaeforthescalingfunctionandwavelet canbeexpressedintermsofthemulti-scalerelationsM1andM2: x = L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k c k x )]TJ/F11 9.9626 Tf 9.962 0 Td [(k ; .7 x = L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k d k x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k : .8 Toobtainaninitialrepresentationofthescalingfunctionandwaveletweonlyneedtodenethe scalingcoecientsckanddk.InthecaseoftheDaubechiesL=4scalingfunctionandwavelet, weexpresstheorthogonalityconditionbyrequiring[14]: d 0 = c 3 ;d 1 = )]TJ/F11 9.9626 Tf 7.749 0 Td [(c 2 ; d 2 = c 1 ;d 3 = )]TJ/F11 9.9626 Tf 7.749 0 Td [(c 0 ; .9 whichissucienttomakethedotproduct x x vanish: x x = h c 0 ;c 1 ;c 2 ;c 3 ih d 0 ;d 1 ;d 2 ;d 3 i = h c 0 ;c 1 ;c 2 ;c 3 ih c 3 ; )]TJ/F11 9.9626 Tf 7.749 0 Td [(c 2 ;c 1 ; )]TJ/F11 9.9626 Tf 7.749 0 Td [(c 0 i = c 0 c 3 )]TJ/F11 9.9626 Tf 9.963 0 Td [(c 1 c 2 + c 2 c 1 )]TJ/F11 9.9626 Tf 9.963 0 Td [(c 3 c 0 =0 : .10 Havingequatedscalingcoecientsforthescalingfunctionandwavelet,wecansolveforthe coecientsbydeningasucientnumberofequationswhichdescribeadesiredbehaviorforour basisfunctions.ThisamountstoaddingconditionsontheregularityofhigherorderDaubechies wavelets[8].Thersttwoconditionswewishtheabovecoecientstosatisfyareaconditiononthe 16 PAGE 21 sizeofthescalingfunction: L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k =0 c k =2 ; .11 thisconditionisactuallyrequiredforanyorthonormalwaveletbasis,seeappendixAandacondition thatshiftedbasisfunctionsaremutuallyorthogonalwheretheysharesupport: h c 0 ;c 1 ;c 2 ;c 3 ; 0 ; 0 ih 0 ; 0 ;c 0 ;c 1 ;c 2 ;c 3 i =0 0 c 0 +0 c 1 + c 0 c 2 + c 1 c 3 +0 c 2 +0 c 3 =0 ; .12 whichisequivalentto: c 0 c 2 + c 1 c 3 =0 : .13 RequirementsfororthogonalityandcoecientsizeholdforallDaubechieswaveletsandscalingfunctions,includingtheHaarwaveletdenedastheDaubechieswaveletwithL=2scalingcoecients. TheadditionalregularityrequirementsfortheDaubechiesL=4waveletbasisarethatwecan expressconstantandlinearfunctionsentirelyintermsofthescalingfunction.Thedotproductof h 1 ; 1 ; 1 ; 1 i withawaveletshouldthenvanish: c 3 )]TJ/F11 9.9626 Tf 9.963 0 Td [(c 2 + c 1 )]TJ/F11 9.9626 Tf 9.962 0 Td [(c 0 =0 ; .14 asshouldthedotproductofthewaveletwithaconstantfunction h 1 ; 2 ; 3 ; 4 i : c 3 )]TJ/F8 9.9626 Tf 9.963 0 Td [(2 c 2 +3 c 1 )]TJ/F8 9.9626 Tf 9.962 0 Td [(4 c 0 =0 : .15 Theabovetwoconditionsallowustoexactlydescribeconstantandlinearfunctionsusingthescaling functionasabasis.HigherorderDaubechieswaveletsareabletoexactlyrepresenthigher-order polynomialsbyrequiringsimilarconstraintsonthescalingcoecients[5].Namely,thescaling 17 PAGE 22 coecientsarerequiredtoobey X k )]TJ/F8 9.9626 Tf 7.749 0 Td [(1 k k m c k =0 .16 forincreasingordersofm. Oncewecalculatevaluesforthesescalingcoecients,wecanthenconstructthescalingfunction andsubsequently,thewavelet.Givenvaluesforthescalingcoecients,onewaytoconstructthe scalingfunctionisviathecascadealgorithm.Thecascadealgorithmbeginswithintegervaluesof thescalingfunctionwhicharegivenexplicitlybythescalingcoecientsandthencalculatesvalues forthescalingfunctionatthehalfintegersusingthemulti-scalerelationsM1andM2. ThecascadealgorithmfortheDaubechiesL=4scalingfunction: [15] Thealgorithmproceedsiteratively,calculatingvaluesateachstepathalfthespacingoftheprevious 18 PAGE 23 step,andthenllinginvaluesforthescalingfunctionusinglinearinterpolationoncetherepresentation issucientlydense. 19 PAGE 24 Chapter3 RepresentingFunctionsand DerivativeOperatorswithaWavelet Basis 3.1ProjectingFunctionsintoaWaveletBasis Givenanexpressionofafunctionasaweightedsuminvolvingawaveletbasis.3,wewantto calculatethecoecients, k .Wecanutilizetheorthogonalityoftheorthonormalwaveletbasis functions, < k j k 0 > = kk 0 ; .1 togenerateaformulaforthe k 'sapplyingthismethodtosineorcosinefunctionsisknownas Fourier'sTrick.Atagivenscale,wecanuse2.3torepresentfx: f x = X k k k x .2 20 PAGE 25 multiplyingbothsidesby l andintegratingyields Z l x f x dx = Z l x X k k k x dx = X k k Z l x k x dx: .3 Applyingtheorthogonalitycondition.1on : X k k Z l x k x dx = X k k l;k = l : .4 Thuswecanndthe k coecientsbyprojectingtheoriginalfunctionintothewaveletbasis: l = Z f x l x dx: .5 Forexample,taking f x =1 givescoecients k = Z 1 k x dx .6 whichisidentically1bythenormalizationcondition: Z x dx =1 : .7 Thus f x =1= X k k k x = X k k x .8 andweseethattoconstruct f x =1 withwaveletswejustaddtogetherwaveletsthataretranslated byk.Thatis,ifweaddtogetherscalingfunctionsastheyaretranslatedbyincrementsk,thenthey resultingfunctionshouldbegintolooklike f x =1 : 21 PAGE 26 Similarly,wecantake f x = x k = Z x k x dx .9 substituting y = x )]TJ/F11 9.9626 Tf 9.962 0 Td [(k = Z y + k y dy = Z y y dy + k Z y dy .10 wherewedenethe p th continuousmoment m p = Z x p x dx: .11 Thus f x = x = X k m 1 + km 0 k x .12 22 PAGE 27 andagainwecanaddtranslatedwaveletswithweights k = m 1 + km 0 toformtheoriginal function: 3.2ContinuousandDiscreteWaveletMoments Evaluatingthebasiscoecientsforawaveletrepresentationofafunctionofteninvolvesprojecting apolynomialinxintoawaveletbasis.Thisrequiresthatweevaluateintegralsoftheform m p = Z x p x dx: .13 wherewecall m p the p th continuousmomentforthescalingfunction x Tobegin,wecancalculatetherstcontinuousmoment, m 1 ,byapplyingthemulti-scalerelation M1 m 1 = Z x x dx = L )]TJ/F7 6.9738 Tf 6.226 0 Td [(1 X k =0 c k Z x x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k dx: .14 23 PAGE 28 Substituting x = z + k = 2 yields m 1 = L )]TJ/F7 6.9738 Tf 6.226 0 Td [(1 X k =0 c k Z z + k 2 z dz 2 = 1 4 L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k =0 c k Z z z dz + k Z z dz ; .15 whichwerecognizeasasumoftherstandzerothcontinuousmoments: m 1 = 1 4 L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k =0 c k [ m 1 + km 0 ] : .16 Nowwehaveanequationforthecontinuousmoment m 1 intermsofsumsoverthecoecients c k Wedenethe p th discretemoment, p as p = 1 2 X k kc k .17 andwerewritetheaboveformulafor m 1 intermsofthesediscretemoments: m 1 = m 1 2 0 + m 0 2 1 : .18 Thenormalizationcondition.7givesus m 0 =1 ,where.11implies 0 =1 ,andwecanrewrite theaboveequationas: m 1 = m 1 2 + 1 2 ; .19 m 1 2 = 1 2 .20 m 1 = 1 .21 andwehaveexpressedanintegralnecessaryforprojectingfunctionsintowaveletspaceasadiscrete sumwhichcanbeevaluatednumerically.Similarresultscanbederivedforanyordercontinuous 24 PAGE 29 moment,withageneralformulagivenby[16]: m p = )]TJ/F8 9.9626 Tf 9.962 0 Td [(2 )]TJ/F10 6.9738 Tf 6.226 0 Td [(p )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 p )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X p 0 =0 p p 0 m p 0 p )]TJ/F10 6.9738 Tf 6.226 0 Td [(p 0 ; .22 wherethe p arecalledthediscretemomentsandaregivenby: p = 1 2 L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k =0 c k k p : .23 3.3ApproximatingProjectionCoecientsusingNumericalQuadrature Whenwelackananalyticdescriptionfortheprojectioncoecientswecanapproximatethemby interpolatingvaluesforthefunctionfxalonganumberofsamplepoints, f x q g q r WecanapproximateanysmoothfunctionasasumofLagrangebasispolynomials.Thatis, f x L r x = r X q =1 f x q l x q .24 wherertheorderoftheinterpolationpolynomialgivesthenumberofLagrangebasispolynomials inthesum.Thebasispolynomialsthemselvesaregivenby: l r;q x = r Y q 0 6 = q x )]TJ/F11 9.9626 Tf 9.963 0 Td [(x q 0 x q )]TJ/F11 9.9626 Tf 9.963 0 Td [(x q 0 : .25 Formonomialsofdegree PAGE 30 Lagrangebasisfunctions: x n = r X q =1 x n q l x q : .26 Theprojectioncoecientsforpolynomialsarealreadygivenbythecontinuouswaveletmoments,but inthecaseofmoregeneralfunctionswewanttobeabletoapproximatethecoecientsusingthe interpolatedvaluesforfxacrossthesamplepoints, f x q g : Z x f x r X q =1 q f x q : .27 Wecancalculatethevalues[17]forthe q atresolutionlevelJas: q = r )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X p =0 2 )]TJ/F9 4.9813 Tf 5.396 0 Td [(J 2 )]TJ/F10 6.9738 Tf 6.226 0 Td [(J p m p p l p r;q k 2 J : .28 Withtheaboveformulaewecanusewaveletstorepresentarbitrarysquare-integrablefunctions, althoughtheirimplementationindescribingpolynomialfunctionsasintheprevioustwosections remainsmuchsimpler. 3.4Finite-DierenceDerivativeOperators WhilethemethodsofthissectionareusedonlyintheFinite-DierenceTime-Domainsimulation ofelectromagneticwavepropagationandnotinMulti-ScaleTime-Domainsimulations,nitedierencemethodsoerasimpleintroductiontoapplyingderivativeoperatorstofunctionsina discretizedspace.Specically,thenite-dierencemethodappliesabandedderivativematrixtoa vectorrepresentingdiscretizedvaluesofsomefunction,fx.Thederivativematrixactsoneach pointinthevectorinsuchawaythatthenewpointmapstotheaccompanyingvalueoff'x.While themethodsofthissectionapplytoasimplediscretizedsamplingofsomefunction,themethod outlinedinsection3.5appliestofunctionsprojectedintoawaveletbasis,changingtheformofthe 26 PAGE 31 associatedderivativeoperators. Wecanapproximatethederivativeofafunctionbycalculatingtheslopeofthesecantlinewhich passesthroughthevaluesofthefunctionattwopointsseparatedbyasmalldistance, x : ForwardDierenceApproximation: f 0 x f x + x )]TJ/F11 9.9626 Tf 9.963 0 Td [(f x x ; .29 BackwardDierenceApproximation: f 0 x f x )]TJ/F11 9.9626 Tf 9.963 0 Td [(f x )]TJ/F8 9.9626 Tf 9.963 0 Td [( x x : .30 Wherethecontinuousrst-orderderivativeoffxdiersfromtheaboveapproximationsonlyinthat x 0 andthesecantlineconvergestoalinetangenttofx. Both.29and.30arerst-orderapproximations{theextenttowhichtheyaccuratelyapproximatethecontinuousderivativef'xvarieslinearlywithspacing x .Toseethis,wecanconsider theTaylorexpansionof f x + x : f x + x = f x + x f 0 x + x 2 2 f 00 x + :::; .31 fromwhichwesubtractfxandre-writethehigher-ordertermsasapolynomialin x : f x + x )]TJ/F11 9.9626 Tf 9.962 0 Td [(f x = x f 0 x + O x 2 : .32 Byrearrangingtheaboveterms,wecanwriteanerrorterm, O x ,whichdescribesthedierence betweenf'xanditsnite-dierenceapproximation.Fortheforwarddierenceapproximationwe have: f 0 x = f x + x )]TJ/F11 9.9626 Tf 9.963 0 Td [(f x x + O x : .33 27 PAGE 32 Higher-orderderivativescanalsobewrittenasnite-dierenceapproximations.Forinstance, wecanapproximatef"xasthedierencebetweentheforward-dierenceandbackward-dierence approximationsoff'x: f 00 x f x + x )]TJ/F10 6.9738 Tf 6.227 0 Td [(f x x )]TJ/F10 6.9738 Tf 11.158 4.832 Td [(f x )]TJ/F10 6.9738 Tf 6.226 0 Td [(f x )]TJ/F7 6.9738 Tf 6.227 0 Td [( x x x = f x + x )]TJ/F8 9.9626 Tf 9.962 0 Td [(2 f x + f x )]TJ/F8 9.9626 Tf 9.962 0 Td [( x x 2 .34 Finite-dierenceapproximationsliketheonesabovehaveadistinctadvantageincomputational simulationssincewecanusethemtoexpressthederivativesofafunctionusingbasicoperations addition,subtraction,multiplication,anddivision.Implementingtheseapproximatederivativeoperatorsisthusrelativelystraightforward.Thenextstepinnumericallyevaluatingagivendierential equationisthentoapplythesederivativeoperatorstoaspatially-discretizedsamplingofagiven function.Thatis,bysamplingafunctionfxatintervalsseparatedby x ,wecanexpressthe derivativeoperatorasamatrixwhichappliesthenite-dierenceapproximationtoeachsamplepoint offx.Discretizingthesecond-orderderivativeoperatorcanbeexpressedasamatrixoftheform: d 2 dx 2 1 x 2 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1000 ::: 0 1 )]TJ/F8 9.9626 Tf 7.749 0 Td [(210 0 01 )]TJ/F8 9.9626 Tf 7.749 0 Td [(21 0 001 )]TJ/F8 9.9626 Tf 7.749 0 Td [(2 0 . . . . 0000 ::: 1 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 .35 which,whenappliedtoanappropriatelysampledfunction,yieldsanumericalapproximationforthe 28 PAGE 33 second-orderderivativeofthefunction: d 2 f dx 2 1 x 2 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1000 ::: 0 1 )]TJ/F8 9.9626 Tf 7.749 0 Td [(210 0 01 )]TJ/F8 9.9626 Tf 7.749 0 Td [(21 0 001 )]TJ/F8 9.9626 Tf 7.749 0 Td [(2 0 . . . . 0000 ::: 1 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 f 1 f 2 f 3 f 4 f n 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 f 1 x 2 f 1 )]TJ/F7 6.9738 Tf 6.227 0 Td [(2 f 2 + f 3 x 2 f 2 )]TJ/F7 6.9738 Tf 6.227 0 Td [(2 f 3 + f 4 x 2 f 3 )]TJ/F7 6.9738 Tf 6.227 0 Td [(2 f 4 + f 5 x 2 f n x 2 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 .36 Tosolveadierentialequation,weneedonlyapplytheappropriatecoecientsandderivativeoperatorstoadiscretizedsamplingofagivenfunction.Forexample,consider )]TJ/F10 6.9738 Tf 8.021 -4.147 Td [(d dx 2 f = f x : 1 x 2 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1000 ::: 0 1 )]TJ/F8 9.9626 Tf 7.749 0 Td [(210 0 01 )]TJ/F8 9.9626 Tf 7.748 0 Td [(21 0 001 )]TJ/F8 9.9626 Tf 7.749 0 Td [(2 0 . . . . 0000 ::: 1 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 f 1 f 2 f 3 f 4 f n 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 f 1 f 2 f 3 f 4 f n 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 : .37 Byexpressingderivativeoperatorsasmatriceswhichactondiscretizedsamplingsofafunction,we converttheproblemofsolvingadierentialequationintoaproblemofnumericallyevaluatinga systemofalgebraicequationswhichisamuchbetter-suitedproblemforcomputationalevaluation. Thisisalsothegoalinprojectingderivativeoperatorintoawaveletbasisasinthefollowingsection, althoughthemachineryismorecomplicated. 3.5ProjectingDerivativeOperatorsintoWaveletSpace Inordertosimulatethetime-dependentbehaviorofasystemitisnecessarytooperateonaspatiallydiscretizedfunctiondescribedasaseriesofcoecientsinawaveletbasisusingsomederivative 29 PAGE 34 operator.Numericallyevaluatingthederivativeofafunctioninvolvesspatiallydiscretizingthederivativeoperator,whichweaccomplishbyprojectingtheoperatorinawaveletbasis.Todothis,we calculatethequantity < j ^ H j k > ,whichinvolvesevaluatingintegralsoftheform r qk = Z x d dx q x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k dx .38 Usingthemulti-scalerelationM1,thiscanbere-expressedasadouble-sumspecifyingthecoecient inthe q th rowand k th columnofaderivative-matrixr: r qk = L )]TJ/F7 6.9738 Tf 6.226 0 Td [(1 X k 0 =0 L )]TJ/F7 6.9738 Tf 6.226 0 Td [(1 X k =0 c k 0 c k Z x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k 0 d dx q x )]TJ/F8 9.9626 Tf 9.963 0 Td [(2 k )]TJ/F11 9.9626 Tf 9.963 0 Td [(k dx: .39 Substituting y =2 x )]TJ/F11 9.9626 Tf 9.962 0 Td [(k 0 yields r qk = L )]TJ/F7 6.9738 Tf 6.226 0 Td [(1 X k 0 =0 L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k =0 c k 0 c k Z y q d dy q y )]TJ/F8 9.9626 Tf 9.963 0 Td [(2 k )]TJ/F11 9.9626 Tf 9.963 0 Td [(k + k 0 1 2 dy .40 resultinginanequationwhichrelatesthecoecientstoeachother: r qk =2 q )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k 0 =0 L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X k =0 c k 0 c k r q; 2 k + k )]TJ/F10 6.9738 Tf 6.226 0 Td [(k 0 : .41 Wethenx r qk byrequiring[18]: q != X k x qk r qk .42 where x qk isamatrixwhose q th rowgivesprojectioncoecientsalongthekcolumnsforthe q th ordermonomialinx: x qk = Z x q k x .43 30 PAGE 35 Chapter4 Maxwell'sEquationsina Schrodinger-likeForm 4.1Maxwell'sEquations TosimulatethetimedependenceofanelectromagneticeldwelooktoMaxwell'sequations,which givepropagationthroughfreespaceoftheelectriceld E andmagneticinduction B as: E = c r B .1 and B = )]TJ/F11 9.9626 Tf 7.749 0 Td [(c r E : .2 Inordertosimulatepropagationthroughadielectricmedium,wemustalsoconsiderhowthematerial respondstothepresenceoftheelectromagneticwavesandaectstheirpropagation.Themeasure ofhowanelectromagneticeldinducesdipolemomentsinthematerialthroughwhichittravelsis 31 PAGE 36 givenbythepolarizationdensity, P : P + P = p 2 E ; .3 wheretheattenuationconstant, ,isameasureofhowtheelddiesoasitpassesthroughthe materialandtheplasmafrequency, p ,describestheoscillationoftheelectricchargedensityin metals.Chapter6describesthesimulationofanelectromagneticwavethroughsolidgold,where thesequantitiesaregiven[19]as p =2.16x 10 15 Hzand =1.75x 10 13 Hz. 4.2ConstructingaHamiltonianOperator Inordertoquicklysimulatethepropagationofanelectromagneticwavewewanttochangethe dierentialequationsdescribingtheeldbehaviorintoasystemofalgebraicequations.Inorderto dothis,weneedtoexpressthetimedependenceoftheelectromagneticeldsintheformofan equationrelatingthetime-andspace-derivativesofsomeeigenvector, : = 2 6 6 6 6 6 6 4 E B Q 3 7 7 7 7 7 7 5 ; .4 wherewedenethevector Q as: Q = 1 p P : .5 Wecanmatchequations.1,.2,.3,and.5todescribethetime-dependenceoftheelds in : E = c r B )]TJ/F11 9.9626 Tf 9.962 0 Td [(! p Q B = )]TJ/F11 9.9626 Tf 7.749 0 Td [(c r E Q = p E )]TJ/F11 9.9626 Tf 9.963 0 Td [( Q 32 PAGE 37 FromthiswecaneasilyconstructaHamiltonianoperator[20]toacton : ^ H = 2 6 6 6 6 6 6 4 0 ic r)]TJ/F11 9.9626 Tf 33.763 0 Td [(i! p )]TJ/F11 9.9626 Tf 7.748 0 Td [(ic r 00 i! p 0 )]TJ/F11 9.9626 Tf 7.749 0 Td [(i 3 7 7 7 7 7 7 5 .6 TheHamiltonianoperatorexpressesthespatialderivativesof inawaywhichallowsustowrite Maxwell'sequationsinaformresemblingthetime-dependentSchrodingerEquation: i d d t = ^ H : .7 Notethatif iszerothentheaboveHamiltonianisHermitianthatis, H y =H,where H y isthe conjugatetransposeofthematrixH.Non-HermitianHamiltoniansappearinthecaseoflossymedia andrequireslightlymorecomplicatedalgorithmstosimulateelectromagneticwavepropagation,as describedinsection5.3. Thegeneralsolutionfor attimetisgivenby: t = e )]TJ/F10 6.9738 Tf 6.227 0 Td [(i t )]TJ/F10 6.9738 Tf 6.227 0 Td [(t 0 ^ H t 0 .8 andsowe'vetranslatedtheproblemofndingthevalueoftheelectric,magnetic,andpolarization eldsatagiventimeintoaproblemofnumericallyevaluatingtheabovematrixexponentialand applyingittosomeinitialvalueof 33 PAGE 38 Chapter5 IterativeAlgorithms Oncewehaveestablishedourwaveletbasis,weborrowfromthematrixmethodsofquantummechanicstoformulateanalgorithmwhichallowsustopropagatestationaryandnon-stationarywaves. Thematrixdevelopedforaparticularprobleminourcase,theHamiltonianderivedfromMaxwell's equationsrepresentstheunderlyingphysicallawswewishtosimulate.TypicallytheHamiltonian matricesusedinevenone-dimensionalproblemsareverylarge;theygrowevenlargerasweaddressadditionalspatialdimensions.SincetheHamiltoniandescribestheoverallpropertiesofour system,andsinceourcalculationswillinvolveapplyingthematrixtoavectorrepresentationofthe electromagneticeldmanytimes,ithelpstoreducethematrixtoaworkableform. Inmostcases,thesizeofthematrixmakesdirectdiagonalizationcomputationallyintractable, soweusealgorithmsdesignedtogiveusasimpliedsquarematrixofdimensionmasopposedto thedimensionoftheoriginalsquareHamiltonianmatrix,n,wherem << n.Thissimpliedmatrix shouldhavesimilarpropertiesasthatofourHamiltonianmatrixi.e.,thedominanteigenvaluesof theHamiltonianmatrixshouldcorrespondtotheeigenvaluesofthesimpliedmatrix.Byreducing thedimensionofthematrixandapproximatingthebehavioroftheoriginalHamiltonian,iterative algorithmstendtoreducetheproblemofdeterminingtheeigenvaluesoftheHamiltonianmatrixfrom taking O n 3 stepstoaproblemrequiring O n 2 steps{asignicantimprovement[21]. 34 PAGE 39 Section5.1demonstrateshowasimpleiterativealgorithmaddressesthegeneralissueofapproximatingtheeigenvectorassociatedtothedominanteigenvalueoftheHamiltonianmatrix,whilethe latersectionsdescribemorecomplicatedalgorithmswhichwecanusetosimulateelectromagnetic wavepropagation. 5.1PowerIteration Oneofthesimplestwaystodeterminethebehavioroflargesparsematricesistoexaminejusttheir dominanteigenvalue.Wecanaccomplishthisusinganiterationtechniqueknownaspoweriteration. Thebasicideaisthatmultiplyingamatrixwitharandomlygeneratedvectoroutputsadierent vectorwithvaluesclosertothatofaneigenvectorcorrespondingtothedominanteigenvalueofthe matrix.Thatis,byiterativelymultiplyingtheoutputofamatrix-vectorproductwiththeoriginal matrixagainandnormalizingtheresultingvectorateachiteration,wecansolveforthedominant eigenvalueofasystem. Wecandemonstratepoweriterationbyndingtheeigenvectorwhichcorrespondstothedominant eigenvalueforanexample3x3matrix.Givensomematrix A = 2 6 6 6 6 6 6 4 231 46 )]TJ/F8 9.9626 Tf 7.748 0 Td [(4 35 )]TJ/F8 9.9626 Tf 7.748 0 Td [(2 3 7 7 7 7 7 7 5 .1 wecandirectlydetermineitseigenvaluesbysolvingthecharacteristicequation: 2 )]TJ/F11 9.9626 Tf 9.962 0 Td [( 31 46 )]TJ/F11 9.9626 Tf 9.962 0 Td [( )]TJ/F8 9.9626 Tf 7.749 0 Td [(4 35 )]TJ/F8 9.9626 Tf 7.749 0 Td [(2 )]TJ/F11 9.9626 Tf 9.963 0 Td [( = 3 )]TJ/F8 9.9626 Tf 9.962 0 Td [(6 2 + )]TJ/F8 9.9626 Tf 9.963 0 Td [(6=0 .2 35 PAGE 40 whichyieldseigenvalues f 1 ; 2 ; 3 g = f 6 ;i; )]TJ/F11 9.9626 Tf 7.748 0 Td [(i g .3 withcorrespondingeigenvectors ~v 1 = 0 B B B B B B @ 1 1 1 1 C C C C C C A ; ~v 2 = 0 B B B B B B @ 65 )]TJ/F8 9.9626 Tf 7.749 0 Td [(44+12 i 2+29 i 1 C C C C C C A ; ~v 3 = 0 B B B B B B @ 65 )]TJ/F8 9.9626 Tf 7.749 0 Td [(44 )]TJ/F8 9.9626 Tf 9.963 0 Td [(12 i 2 )]TJ/F8 9.9626 Tf 9.963 0 Td [(29 i 1 C C C C C C A : .4 Thepowermethodshouldtakeasinputanyvectorin C n andoutputavectorwhichapproximates thevectorassociatedwiththedominanteigenvalueoftheoriginalmatrix.Thatis,since k 1 k > k 2 k ; k 3 k ,weexpectthatpoweriterationwillresultinprogressivelybetterapproximationsofthe eigenvector ~v 1 .Toseeifthepoweriterationmethodyieldssimilarresults,wecouldchoosean arbitraryvector, ~u ,andrepeatedlymultiplyitbythematrixA: A 2 u = 2 6 6 6 6 6 6 4 1929 )]TJ/F8 9.9626 Tf 7.749 0 Td [(12 2028 )]TJ/F8 9.9626 Tf 7.749 0 Td [(12 2029 )]TJ/F8 9.9626 Tf 7.749 0 Td [(13 3 7 7 7 7 7 7 5 2 6 6 6 6 6 6 4 7 33 )]TJ/F8 9.9626 Tf 7.749 0 Td [(6 3 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 4 1162 1136 1175 3 7 7 7 7 7 7 5 ; .5 A 4 u = 2 6 6 6 6 6 6 4 7011015 )]TJ/F8 9.9626 Tf 7.749 0 Td [(420 7001016 )]TJ/F8 9.9626 Tf 7.749 0 Td [(420 7001015 )]TJ/F8 9.9626 Tf 7.749 0 Td [(419 3 7 7 7 7 7 7 5 2 6 6 6 6 6 6 4 7 33 )]TJ/F8 9.9626 Tf 7.749 0 Td [(6 3 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 4 40922 40948 40909 3 7 7 7 7 7 7 5 ; .6 A 8 u = 2 6 6 6 6 6 6 4 9079011316455 )]TJ/F8 9.9626 Tf 7.749 0 Td [(544740 9079001316456 )]TJ/F8 9.9626 Tf 7.749 0 Td [(544740 9079001316455 )]TJ/F8 9.9626 Tf 7.749 0 Td [(544739 3 7 7 7 7 7 7 5 2 6 6 6 6 6 6 4 7 33 )]TJ/F8 9.9626 Tf 7.749 0 Td [(6 3 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 4 53066762 53066788 53066749 3 7 7 7 7 7 7 5 : .7 Intheaboveexampleitisalsoeasytoseethatthecolumnsoftheresultingmatricesconvergeto 36 PAGE 41 theeigenvectorcorrespondingtothedominanteigenvalueoftheoriginalmatrix. Notethatalthoughthedirectcomputationoftheeigenvaluesfortheabove3x3matrixissimple inthiscase,theoperationsinvolvedbecomemuchmoredemandingasweincreasethedimensionsof thematrix.Itisalsonotalwaysthecasethatthepowermethodwillconvergeasquicklytothedesired eigenvectorasitdidintheaboveexample.Althoughtheyaremorecomplicatedtoimplementand mayrequiremorememory,othermethodsincludingKrylovmethods,asdescribedbelowimprove uponthebasictechniquesetforthbypoweriterationbyconvergingfastertothedesiredeigenvectors andallowingustoapproximateeigenvaluesbesidesthesinglelargesteigenvalue. 5.2KrylovSubspaceTechniques Thecentralproblemwhichwewishtosimplifyisthatofcomputingthematrixexponentialquantities whichdescribethestateofsolutionstoMaxwell'sequationsatagiventime.Onewaytosolvefor t isbyexpanding.8inaTaylorseries: t = e )]TJ/F10 6.9738 Tf 6.227 0 Td [(i t )]TJ/F10 6.9738 Tf 6.227 0 Td [(t 0 ^ H t 0 = 1 X n =0 )]TJ/F11 9.9626 Tf 7.748 0 Td [(i t n n ^ H n t 0 .8 wherethechiefcomputationaltaskbecomesevaluatingthegivenmatrix-vectorpairs.Whilewecan approximatesolutionsto t bychoosingasuitablysmalltimespacing, t ,andtruncatingthe seriesafteranitenumberofterms[22],theHamiltonianmatrixinvolvedtendstobelargeand sparse,whichmakesevaluatingthematrix-vectorproductscumbersome. Alternately,wecansignicantlyreducethecomputationalcomplexityofnding Psi t byapproximatingtheHamiltonianoperatorandapplyingittoasuitablebasisvector[1].Weapproximatethe Hamiltonianbytakingitsprojectioninthesubspace K m whichspanstherstm-1matrix-vector productsgeneratedbypoweriteration: K m = span f 0 ; ^ H 0 ; ^ H 2 0 ;:::; ^ H m )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 0 g : .9 37 PAGE 42 ThisisthemainideabehinditerativealgorithmsknownasKrylovmethods.AlthoughsomeprecautionsmustbetakentopreventthebasisoftheabovesubspacecalledtheKrylovsubspacefrom becominglinearlydependent,wecanusethevectorsfromthesubspacetoformamatrixwhose eigenvaluesapproximatethedominanteigenvaluesoftheoriginalHamiltonian. 5.3TheLanczosAlgorithm Wecanorthonormalizeeachnewvectorinthesubspaceagainstitspredecessorsusingaprocess similartotheGram-Schmidtprocedure.TheLanczosalgorithmusestheseorthonormalizedvectors toproduceatridiagonalmatrixwithdimensionequaltothatoftheKrylovsubspace,theeigenvalues ofwhichapproximatethedominanteigenvaluesoftheoriginalHamiltonianmatrix.Fornon-Hermitian matriceswecanapplytheBi-orthogonalLanczosalgorithm[23],therststepofwhichconstructsa Krylovsubspacefromaninitialpairofnormalizedvectors v 1 w 1 : k ~v 1 k = k ~w 1 k =1 : .10 Whensimulatingthepropagationofelectromagneticwaves,theseinitialvectorsrepresenttheinitial conditionoftheelectromagneticeld.ThebiorthogonalLanczosalgorithmthendenessubsequent valuesof v j and w j bytakingmatrix-vectorproductsoftheHamiltonianmatrixandthepreceding vectors: ~ v j +1 = ^ H ~v j )]TJ/F11 9.9626 Tf 9.962 0 Td [( j ~v j )]TJ/F11 9.9626 Tf 9.963 0 Td [( j ~v j )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 .11a ~ w j +1 = ^ H y ~w j )]TJ/F11 9.9626 Tf 9.963 0 Td [( j ~w j )]TJ/F11 9.9626 Tf 9.962 0 Td [( j ~w j )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 .11b 38 PAGE 43 wheretheconstants j j ,and j aregivenby: j = ~w y j ^ H ~v j .12a j = k ~w j k .12b j = k ~v j k ; .12c andwheretheresultingvectorsarenormalizedbeforethenextiteration: ~v j +1 = ~ v j +1 n +1 .13a ~w j +1 = ~ w j +1 n +1 : .13b OncetheKrylovsubspaceislargeenough,thealgorithmterminatesandweconstructatridiagonal matrix, T m ,whoseeigenvaluesapproximatethedominanteigenvaluesoftheoriginalHamiltonian matrix: T m = 0 B B B B B B B B B B B B B B @ 1 2 0 ::: 0 2 2 3 . 0 3 . 0 . . . m 0 ::: 0 m m 1 C C C C C C C C C C C C C C A : .14 Analnote:wecanformtheKrylovsubspacefromasmanymatrix-vectorproductsaswelike uptothenumberofthedimensionoftheoriginalHamiltonian,atwhichpointtheeigenvalues oftheresultingmatrixexactlyequaltheeigenvaluesoftheHamiltonian.However,inpractice theeigenvaluesofthetridiagonalmatrixconvergequicklytothoseoftheHamiltonianandsothe dimension,m,oftheunderlyingKrylovsubspaceisusuallyaround10. 39 PAGE 44 Chapter6 ExampleSimulation Tosimulatethepropagationofanelectromagneticeldinone-dimensionwedeneasourcefunction fortheelectriceld: whichwediscretizebycalculatingtheprojectioncoecientsofthewaveletbasisusedtorepresent thefunction.Theresultisaspatiallydiscretizedelectriceldvector: 40 PAGE 45 . ThenextstepistogenerateaKrylovbasisbyapplyingtheHamiltonianmatrix.6totheabove eldvector.WethenruntheLanczosalgorithm,usingthevectorsintheKrylovsubspacetogenerate thecoecientsofatridiagonalmatrixwhichismuchsmallerthantheoriginalHamiltonianmatrix: 41 PAGE 46 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ 0 28.9456794 i 0 : 501295700000000 0 : 5012957 )]TJ/F8 9.9626 Tf 7.749 0 Td [(0 : 0172822 i 0 : 49824050000000 00 : 4982405 )]TJ/F8 9.9626 Tf 7.748 0 Td [(0 : 0002069 i 0 : 0090538000000 000 : 00905380 : 0185326 i 0 : 716464800000 0000 : 7164648 )]TJ/F8 9.9626 Tf 7.749 0 Td [(0 : 0354838 i 0 : 00118500000 00000 : 001185010 : 5274038 i 14 : 6534813000 0000014 : 6534813 )]TJ/F8 9.9626 Tf 7.749 0 Td [(15 : 6032195 i 4 : 9244249 i 00 0000004 : 9244249 i 5 : 0589165 i 0 : 0043963 i 0 00000000 : 0043963 i 0 : 0747066 i 0 : 7685698 000000000 : 7685698 )]TJ/F8 9.9626 Tf 7.748 0 Td [(0 : 0893247 i 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A Thetridiagonalmatrixofdimensionm=10generatedbytheLanczosalgorithmfortheHamiltonianoperatordenedin.6.Thismatrixisactually forpropagationthroughlosslessmedia,ittsonapagebetterthantheoneusedforthelossymediasimulation.Thedimensionoftheoriginal discretizedHamiltonianoperatorcanbemanyordersofmagnitudelarger,resultingincumbersomecomputationalstepsin.8. 42 PAGE 47 WecanthensimplyusethismatrixinplaceoftheoriginalHamiltoniantocalculatethematrixexponentialquantities.8andestimatethevaluesoftheelectromagneticeldsatagiventime. Foralossygolddielectricfunction,wecandemonstratethedropointhemagnitudeofthe electromagneticelds: Thedecreaseinthemagnitudeoftheelectromagneticsupervector, ,iscalculatedaboveatdiscrete timesteps t.Wecanalsocalculatevaluesfortheindividualeldsateachofthesetimesteps: 43 PAGE 48 44 PAGE 49 Chapter7 Conclusion Thewaveletbasesdescribedinchapters2and3haveavarietyofapplications,includingsignal processingandmodellingpartialdierentialequations.Theiterativealgorithmsdescribedinchapter 5areembeddedinthesolutionoftime-steppingpropagationmethodsforelectromagneticwavesas describedinchapter4.Specically,byapplyingtheBiorthogonalLanczosAlgorithmtothenonHermitianHamiltoniandescribingalossydielectricfunction,itispossibletoeectivelysimulatethe propagationanddecreaseinmagnitudeoftheelectric,magneticandpolarizationeldsthrougha lossymedium.Inall,theMulti-ScaleTime-Domainmethoddevelopedallowsforthesimulationof electromagneticwavesthroughlossymediaandcanapplytomorecomplicatedsimulations,including hopefullythesimulationofelectromagneticwavespropagatingacrossadiscontinuousdielectric functionorinthepresenceofcurrentandchargesources-scenariostowhichtheFinite-Dierence Time-Domainispoorlysuitedduetothedicultiesinvolvedinmodellingthephysicalpropertiesof thesesystemsusingFourierbasisfunctions. 45 PAGE 50 AppendixA SizeRequirementonScaling Coecients Combiningtheorthogonalityrequirement.1ofanorthonormalwaveletbasiswiththemulti-scale relationM1allowsustoderiveaconstraintonthecombinedsizeofthescalingcoecients.By substitutingy=x+k,wecanchangetheformoftheorthogonalitycondition: < k j k 0 > = Z x )]TJ/F11 9.9626 Tf 9.962 0 Td [(k x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k 0 dx = Z y y + k )]TJ/F11 9.9626 Tf 9.963 0 Td [(k 0 dy = kk 0 ; A.1 whereletting k 0 =0 yields: Z x x + k dx = k; 0 A.2 46 PAGE 51 Byscalingthebasisfromj=0toj=1weintroduceanequationinvolvingthecoecients, c l and c m : L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X l =0 c l L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X m =0 c m Z x )]TJ/F11 9.9626 Tf 9.962 0 Td [(l [ x )]TJ/F11 9.9626 Tf 9.963 0 Td [(k ] )]TJ/F11 9.9626 Tf 9.963 0 Td [(m dx = k; 0 A.3 Substitutingy=2x-l,weobtain: L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X l =0 c l L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X m =0 c m Z y y )]TJ/F8 9.9626 Tf 9.963 0 Td [(2 k )]TJ/F11 9.9626 Tf 9.963 0 Td [(m + l dy 2 = k; 0 A.4 Imposingtheorthogonalityconditionontheintegralagainyields: Z y y )]TJ/F8 9.9626 Tf 9.962 0 Td [(2 k )]TJ/F11 9.9626 Tf 9.962 0 Td [(m + l dy = m +2 k;l ; A.5 andthus k; 0 = 1 2 L )]TJ/F7 6.9738 Tf 6.226 0 Td [(1 X l =0 c l L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X m =0 c m m +2 k;l = 1 2 L )]TJ/F7 6.9738 Tf 6.227 0 Td [(1 X m =0 c m +2 k c m A.6 Thisyieldsaconditiononthewaveletcoecients: L )]TJ/F7 6.9738 Tf 6.226 0 Td [(1 X m =0 c m +2 k c m =2 k; 0 : A.7 Or,letting k =0 L )]TJ/F7 6.9738 Tf 6.226 0 Td [(1 X m =0 c 2 m =2 ; A.8 whichimplies.11. 47 PAGE 52 Bibliography [1]KurtBusch,JensNiegemann,MartinPototschnig,andLashaTkeshelashvili.Akrylov-subspace basedsolverforthelinearandnonlinearmaxwellequations. pss-b ,244:3479{3496,September2007. [2]AllenTaoveandSusanC.Hagness. ComputationalElectrodynamics:TheFinite-Dierence Time-DomainMethod .ArtechHouse,secondedition,2000. [3]RobiPolikar.Cosinefouriertransform,1994.CosineFourierTransformimage. [4]Omegatron.Haarwavelet,June2006.Haarwaveletimage. [5]GilbertStrang.Waveletsanddilationequations:Abriefintroduction. SiamReview ,31:613{ 627,1989. [6]RoyHaandJustinRomberg.Thehaarscalingfunctionandwaveletsatdierentresolutions, July2006.Haarscalingfunctionandwaveletsimage. [7]GeorgeOppenheim. WaveletsandtheirApplications .ISTELtd,2007. [8]M.PollicottandH.Weiss.Howsmoothisyourwavelet?waveletregularityviathermodynamic formalism.June2005. [9]AlfredHaar.Zurtheoriederorthogonalenfunktionensysteme. Math.Ann. ,69:331{371,1910. [10]LutzL.Daubechies4tapwaveletandscalingfunctions,June2005.DaubechiesL=4image. 48 PAGE 53 [11]LutzL.Daubechies12tapwaveletandscalingfunctions,June2005.DaubechiesL=12image. [12]LutzL.Daubechies20tapwaveletandscalingfunctions,June2005.DaubechiesL=20image. [13]FilipWasilewski.Symlet6waveletandscalingfunctions,2008.Symmletscalingfunctionand waveletimage. [14]GilbertStrang.Wavelets. AmericanScientist ,82:250{255,April1994. [15]PhilSchniter.Steps1-5,10ofthecascadealgorithm,June2009.CascadeAlgorithmimage. [16]BruceR.Johnson,JasonP.Modisette,PeterJ.Nordlander,andJamesL.Kinsey.Quadrature integrationfororthogonalwaveletsystems. JournalofChemicalPhysics ,110:8309{8317, May1999. [17]BruceR.JohnsonandJamesL.Kinsey.Quadraturepreltersforthediscretewavelettransform. IEEETransactionsonSignalProcessing ,48:873{875,March2000. [18]G.Beylkin.Ontherepresentationofoperatorsinbasesofcompactlysupportedwavelets. SIAM JournalonNumericalAnalysis ,29:1716{1740,December1992. [19]RamiroAcevedo,RichardLombardini,MatthewA.Turner,JamesL.Kinsey,andBruceR. Johnson.Quantumandelectromagneticpropagationwiththeconjugatesymmetriclanczos method. TheJournalofChemicalPhysics ,128,2008. [20]AndreiG.BorisovandSergeiV.Shabanov.Lanczospseudospectralmethodforinitial-valueproblemsinelectrodynamicsanditsapplicationstoioniccrystalgratings. JournalofComputational Physics ,209:643{644,2005. [21]GeorgeW.CollinsII. FundamentalNumericalMethodsandDataAnalysis .2003. [22]CarlP.Goodrich,RichardLombardini,RamiroAcevedo,andBruceR.Johnson.Non-hermitian propagationtechniquesinamulti-resolutionwaveletbasis.2008. [23]RolandW.Freund.Lanczos-typealgorithmsforstucturednon-hermitianeigenvalueproblems. 49 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||