时域有限差分法的MATLAB应用 值得拥有！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！Contents1 Numeric artifacts1.1 Introduction1.2 Finite Precision81.3 Symbolic manipulation2 Brief Review of Electromagnetics132.1 Introduction2.2 Coulomb,s law and electric Ficld2.3 Electric Flux Density2. 4 Static electric fields172.5 Gradient Divergence and Cur182.6 Laplacian212. 7 Gauss's and Stokes' Theorems242.8 Electric Field boundary conditions252.9 Conductivity and perfect Electric Conductors252.10 Magnetic fields262.11 Magnetic Field boundary conditions272.12 Summary of static Fields272. 13 Time varying fields282. 14 Summary of Time-Varying Fields292. 15 Wave Equation in a Source-Free Region292.16 One-Dimensional Solutions to the Wave Equation303 Introduction to the fdtd method333.1 Introduction333.2 The Yee algorithm343. 3 Update equations in ID353.4 Computcr Implementation of a Onc-DimcnsionalFDTD SImulation393.5 Bare-Bones simulation413.6 PMC Boundary in One Dimension443.7 Snapshots of the Field453. 8 Additive source483.9 Terminating the Grid503. 10 Total-Field/Scattered-Field Boundary53CONTENTS3.11 Inhomogeneities3.12 LosSy Material664 Improving the Fdtd code754.1 Introduction754.2 Arrays and Dynamic Memory Allocation754.3 Macros774.4 Structures804.5 Improvement Number One864.6 Modular design and initialization Functions904.7 Improvement Number Two954.8 Compiling modular Code1024.9 Improvement Number Three..1035 Scaling FDTD Simulations to Any Frequency1155.1 Introduction1155.2 Sources155.2.1 Gaussian pulse1155.2.2 Harmonic Sources1165.2.3 The Ricker Wavelet,,,,,,,,1175.3 Mapping Frequencies to Discrete Fourier Transforms1205.4 Running discrete Fourier Transform (DFT)1215.5 Real Signals and DFTs.,,.,,,,1235.6 Amplitude and Phase from Two Time-Domain Samples1265.7 Conductivity1285. 8 Transmission Coefficient for a Planar Interface........1325.8.1 Transmission through Planar Interface1345.8.2 Mcasuring thc Transmission Cocfficicnt USing FDtD1356 Differential-Equation Based ABCs1456.1 Introduction1456.2 The Advection equation1456.3 Terminating the grid1466.4 Implementation of a First-Order abc1486.5 ABC Expressed Using operator notation.,1536.6 Second-Order Abc56.7 Implementation of a second-Order abc..1587 Dispersion, Impedance, Reflection, and Transmission1617.1 Introduction1617.2 Dispersion in the Continuous World1617.3 Harmonic Representation of the FDTD Method1627. 4 Dispersion in the fdtd grid,,,,1657.5 Numeric Impedance1697.6 Analytic FdtD Reflection and Transmission Coefficients169CONTENTS7.7 Reflection from a Pec1737.8 Interface Aligned with an Electric-Field node1758 Two-Dimensional fDtD Simulationsl818.1 Introduction,18l8.2 Multidimensional arrays1818.3 Two Dimensions: TM2 Polarization1858.4 TM2 Example...1898. 5 The tfSF Boundary for TM Polarization2028.6 TM TFSF Boundary example2088.7 TE Polarization2208. 8 PECS in te and TM2 simulations2248.9 TE Example.2279 Three-Dimensional fdtd2419.1 Introduction2419.2 3D Arrays in C..2419.3 Governing Equations and the 3D Grid244943 D Examplc·2529.5 TFSF Boundary,,,,,,,.2679.6 TFSF Demonstration..2729.7 Unequal Spatial Steps28210 Dispersive material28910.1 Introduction28910.2 Constitutive Relations and Dispersive media29010. 2. Drude materials29110.2.2 Lorentz material29210.2. 3 Debye material29310.3 Debye Materials Using the ADE Method29410.4 Drude Materials Using the ADE Method29610.5 Magnetically Dispersive material29810.6 Piecewise linear recursive convolution30110.7 PLRC for Debye material.30511 Perfectly Matched Layer30711. 1 Introduction30711. 2 Lossy layer, ID30811.3 Lossy Layer, 2D31011. 4 Split-Field Perfectly Matched Layer31211.5 Un-Split Pml...31511. 6 FDTD Implementation of Un-Split PML..318CONTENTS12 Acoustic fdtd simulations32312.1 Introduction32312.2 Governing FDTD Equations32512. 3 Two-Dimensiona11 Impldementation32813 Parallel processing33113. 1 Threads33113.2 Thread Examples33313.3 Message Passing Interface34013.4 Open MPI basics34113.5 Rank and size34313.6 Communicating between Processes34414 Near-to-Far-Field Transformation35114.1 Introduction35114.2 The Equivalence Principldle35114.3 Vector potentials..35214.4 Electric Field in the Far -Field.,,,,,,,,,,,,,,,,,,,,,,.35914.5 Simpsons Composite Intcgration36314.6 Collocating the Electric and Magnetic Fields: The Geometric Mean . .......36314.7 NTFF Transformations Using the Geometric Mean366147.1 Double-Slit radiation36614.7.2 Scattering from a Circular cylinder37014.7.3 Scattering from a Strongly Forward-Scattering Sphere371a Construction of Fourth-Order central differencesA377B Generating a Waterfall Plot and animationB.379C Rendering and animating two-Dimensional DataC383D NotationD387E PostScript PrimerE389E1 IntroductionE.389E. 2 The PostScript File,,,,,,,,E.390E 3 PostScript Basic commandsE.390Index403Chapter 1Numeric artifacts1.1 IntroductionVirtually all solutions to problems in electromagnetics require the use of a computer. Even whenan analytic or " closed form'"solution is available which is nominally exact one typically must usea computer to translate that solution into numeric values for a given set of parameters. Becauseof inherent limitations in the way numbers are stored in computers, some errors will invariably bepresent in the resulting Solution. These crrors will typically be small but they are an artifact aboutwhich one should be aware. Here we will discuss a few of the consequences of finite precisionLater we will be discussing numeric solutions to electromagnetic problems which are basedon the finite-difference time-domain(FDTD) method. The fdtd method makes approximationsthat force the solutions to be approximate, i.e., the method is inherently approximate. The resultsobtained from the FDtd method would be approximate even if we used computers that offeredinfinite numeric precision. The inherent approximations in the FDTD method will be discussed insubsequent chaptersWith numerical methods there is one note of caution which one should always keep in mindProvided the implementation of a solution does not fail catastrophically, a computer is alwayswilling to give you a result. You will probably find there are times when, to get your programsimply to run, the debugging process is incredibly arduous. When your program does run, thenatural assumption is that all the bugs have been fixed. Unfortunately that often is not the caseGetting the program to run is one thing, getting correct results is another. And, in fact, gettingaccurate results is yet another thing-your solution may be correct for the given implementationbut the implementation may not be one which is capable of producing sufficiently accurate resultsTherefore, the more ways you have to test your implementation and your solution, the better. Forexample, a solution may be obtained at one level of discretization and then another solution usinga finer discretization. If the two solutions are not sufficiently close, one has not yet converged tothe"true" solution and a finer discretization must be used or, perhaps, there is some systemic errorin the implementation. The bottom line: just because a computer gives you an answer does notmean that answer is correctLecture notes by John Schneider. numerical-issues tex8CHAPTER 1. NUMERIC ARTIFACTS1.2 Finite PrecisionIf we sum one-eleventh eleven times we know that the result is one, i.e., 1 /11+1/11+1/11+1/11+1/11+1/11+1/11+1/11+1/11+1/11=1. But is that truc on a computcr? Considerthe c program shown in Program 1.1Progran11 onefleventh.c: Test if1/11+1/11+1/11+1/11+1/11+1/11+1/11+1/11+1/11+1/1 1 equals1*Is summing 1./11. ten times --1.0?*,#include 4 int main() t5float ai7a=1.0/11.0;if(a t a+ a t a+ a t a t a+ a t a+ a+ a==l0)printf (" Equal.\n");print f("Not equal.\n")ireturn 0iIn this program the float variable a is set to one-eleventh. In line g the sum of eleven a's iscompared to one. If they are equal, the program prints Equal"but prints Not equal"otherwiseThe output of this program is"Not equal. Thus, to a computer(at least one running a languagetypically used in the solution of electromagnetics problems), the sum of one-eleventh eleven timesis not equal to one. It is worth noting that had line 9 been written a=1/1li, a would have been setto zero since integer math would be used to evaluate the division. By using a =1.0 /110ithe computer uses floating-point mathThe floating-point data types in C or Fortran can only store a finite number of digits. Onmost machines four bytes (32 binary digits or bits) are used for singlc-prccision numbcrs andeight bytes(64 digits)are used for double precision. Returning to the sum of one-elevenths, asan extreme example assumed that a computer can only store two decimal digits. One eleventh isequal to 0.09090909... Thus, to two decimal places one-eleventh would be approximated by 0.09Summing this eleven times yields0.09+0.09+0.09+0.09+0.09+0.09+0.09+0.09+0.09+0.09+0.090.99which is clearly not equal to one. If the number is stored with more digits, the result becomescloser to onc, but it never gcts there. Both the decimal and binary floating-point representationof one-eleventh have an infinite number of digits. Thus, when attempting to store one-eleventh1.2.FⅠ NITE PRECISIONin a computer the number has to be truncated so that the computer stores an approximation ofone-eleventh. Because of this truncation summing one-eleventh eleven times does not yield oneSince 1/10 is equal to 0.1, it might appear this number can be stored with a finite numberof digits. Although one-tenth has a finite number of digits when written in base ten(decimalrepresentation), it has an infinite number of digits when written in base two(binary representation)In a floating-point decimal number each digit represents the number of a particular power often. Letting a blank represent a digit, a decimal number can be thought of in the follow way10110010-110-210-310-4Each digits tells how many of a particular power of 10 there is in a number. The decimal pointserves as the dividing line between negative and non-negative exponents. Binary numbers aresimilar except each digit represents a power or two2The base-ten number 0. 110 is simply 1 x 10-1. To obtain the same value using binary numbers wehave to take 2-4+2-5+2-8+2-9+..., i. e, an infinite number of binary digits. another way ofwriting this is0.10=0.0001100110011001100110011..2As before, when this is stored in a computer, the number has to be truncated. The stored value isno longer precisely equal to one-tenth. Summing ten of these values does not yield one(althoughthe difference is very small)The details of how floating-point values are stored in a computer are not a primary concernHowever, it is helpful to know how bits are allocated Numbers are stored in exponential form andthe standard allocation of bits istotal bits sign ma exponensingle precision328double precision6452Essentially the exponent gives the magnitude of the number while the mantissa gives the digits ofthe number-the mantissa determines the precision. The more digits available for the mantissathe more precisely a number can be represented. Although a double-precision number has twiceas many total bits as a single-precision number, it uses 52 bits for the mantissa whereas a singleprecision number uses 23. Therefore double-precision numbers actually offer more than twice theprecision of single-precision numbers. A mantissa of 23 binary digits corresponds to a little lessthan seven decimal digits. This is because 2-3 is 8, 388, 608, thus 23 binary digits can representnumbers between 0 and 8, 388, 607. On the other hand, a mantissa of 52 binary digits correspondsto a value with between 15 and 16 decimal digits(2 =4, 503, 599, 627, 370, 496)For the exponent, a double-precision number has three more bits than a single-precision num-ber. It may seem as if the double-precision exponent has been short-changed as it does not havetwice as many bits as a single-precision number. However, keep in mind that the exponent represents the size of a number. Each additional bit essentially doubles the number of valucs that can berepresented. If the exponent had nine bits it could represent numbers which were twice as largeCHAPTER 1. NUMERIC ARTIFACTSas single-precision numbers. The three additional bits that a double-precision number possessesallows it to represent exponents which are eight times larger than single-precision numbers. Thistranslates into numbers which are 256 times larger(or smaller)in magnitude than single-precisionnumbersConsider the following equation+6=aFrom mathematics we know this equation can only be satisfied if b is zero. However, using com-puters this equation can be true, i.e., b makes no contribution to a, even when b is non-zeroWhen numbers are added or subtracted, their mantissas are shifted until their exponents areequal. At that point the mantissas can be directly added or subtracted. However, if the differencethe exponents is greater than the length of the mantissa, then the smaller number will not haveany affect when added to or subtracted from the larger number. The code fragment shown inFragment 1.2 illustrates this phenomenonFragment 1.2 Codc fragment to test if a non-zero b can satisfy the cquation a +b=af⊥○ataCia+ bwhile(c ! a) ib=b/2.0;10 printf("12g12g12g\n", a,b, c)Here a is initialized to one while b is set to one-half the variable c holds the sum of a and bThe while-loop starting on line 5 will continue as long as c is not equal to a. In the body of theloop b is divided by 2 and c is again set equal to a +b. If the computer had infinite precision, thiswould be an infinite loop The value of b would become vanishingly small, but it would never bezero and hence a+ b would never equal a. However, the loop does terminate and the output of theprint f() statement in line 10 is5.96046e-08This shows that both a and c are unity while b has a value of 5. 96046x 10. note that this valueof b corresponds to 1 x 2-24. When b has this value, the exponents of a and b differ by more than23(ais1×20)One more example serves to illustrate the less-than-obvious ways in which finite precision cancorrupt a calculation. Assume the variable a is set equal to 2. Taking the square root of a and thensquaring a should yield a result which is close to 2(ideally it would be 2, but since v2 has aninfinite number of digits, some accuracy will be lost). However, what happens if the square root is