Introduction to Partial Differential Equations with Physical Examples, Lecture notes of Differential Equations

Download Introduction to Partial Differential Equations with Physical Examples and more Lecture notes Differential Equations in PDF only on Docsity! Partial differential equations This chapter is an introduction to PDE with physical examples that allow straightforward numerical solution with Mathemat- ica. Methods of solution of PDEs that require more analytical work may be will be considered in subsequent chapters. General facts about PDE Partial differential equations (PDE) are equations for functions of several variables that contain partial derivatives. Typical PDEs are Laplace equation ∆φ@x, y, …D  0 (D is the Laplace operator), Poisson equation (Laplace equation with a source) ∆φ@x, y, …D  f@x, y, …D, wave equation ∂t 2φ@t, x, y, …D − c2 ∆φ@t, x, y, …D  0, heat conduction / diffusion equation ∂tφ@t, x, y, …D − κ∆φ@t, x, y, …D  0, Schrödinger equation  ∂tφ@t, x, y, …D + Ha∆ + bf@x, y, …DL φ@t, x, y, …D  0, etc. There are systems of PDEs and nonlinear PDEs. Solution of PDEs have more freedom than those of ODEs because integration "constants" are in fact functions. For instance, the general solution of the second-order PDE ∂x,yf@x, yD  0 is f@x, yD = F@xD + G@yD, where F@xD and G@yD are arbitrary functions. The solution of the first-order PDE ∂tf@t, xD − v ∂xf@t, xD  0 is f@t, xD = g@x − vtD that describes a front of arbitrary shape moving in the positive direction if v > 0. General analytical solutions of PDEs are available only in the simplest cases, and because of this freedom, they do not yet solve the problem. The actual form of the solution is defined by the symmetry of the problem (if it exists) and boundary conditions. If one of the variables is time, one usually speaks of initial conditions set at the initial time and of the boundary conditions for spatial variables. If there are initial conditions but no final conditions, the problem is evolutionary and it can be solved numerically starting from the initial conditions and increasing time. The most efficient method to do it is the so-called "method of lines" used by Mathematica. First, the problem is discretized in spacial variables and spatial derivatives are approximated by differences. This reduces the PDE to a system of ODEs in time. Then the resulting system of ODEs is solved by one of high-performance ODE solvers. In Mathematica, PDEs, as well as ODEs, are solved by NDSolve. However, currently Mathematica can only solve problems with a rectangular spatial region. For a typical PDE, in these cases one also can find an analytical solution. Mathematically, there is no difference between the time and other variables. If for a particular spacial variable a boundary condition is set only at one end of the interval, one can treat this variable as time and the problem is evolutionary. Mathemat- ica recognizes this situation and finds the solution. If boundary conditions are set at all ends of the interval (or infinity) NDSolve does not find the solution and other methods have to be used. Most time-independent problems are like that. Partial differential equations in physics In physics, PDEs describe continua such as fluids, elastic solids, temperature and concentration distributions, electromag- netic fields, and quantum-mechanical probabilities. Below we review these equations. ü Continuity equation Any conserved quantity described by the density r and the corresponding current j satisfies the continuity equation (1) ∂ρ ∂t + divj  0. This is the continuity equation in the differential form. If „V is an infinitesimal volume, ∑tr„V is the rate of change of the amount (corresponding to r) within this volume, whereas divj „V is the flux out of this volume. If the flux is positive, the amount decreases. The integral form of the continuity equation can be obtained by integrating over a finite volume V ∂ ∂t ‡ ρ V  − ‡ divj V. At the right, the integral over the volume can be expressed via a surface integral, the flux of j, according to the Gauss theorem ∂ ∂t ‡ ρ V  − ‡ j ⋅ S. As there are many types of densities in physics (mass density, energy density, charge density, etc.) there are also many types of continuity equations. ü Fluid dynamics Motion of a fluid satisfies the continuity equation (1), there r is the mass density and the mass current j density is given by j = ρv. Liquids, unlike gases, are practically uncompressible, r = const, thus the continuity equation simplifies to divv  0. Dynamics of a non-viscous fluid is described by the Euler equation ∂v ∂t + Hv ∇L v  − ∇P ρ + g, where P is pressure and g is the gravity acceleration. Minus in front of the pressure gradient shows that the liquid accelerates in the direction of a smaller pressure. The kinematic (streaming) term Hv “L v accounts for the change of the velocity at a given point as a result of the media motion and makes the Euler equation nonlinear. Flow of a viscous fluid obeys the Navier-Stokes equation. For an incompressible fluid it has the form 2 Mathematical_physics-13-Partial_differential_equations.nb ∂e ∂t + div q  0, where e is the density of the internal energy and q ≡ je is its current. If the body is incompressible, the change of the internal energy according to thermodynamics is due to the heat received, „e = dQ. On the other hand, the heat received is related to the temperature change via heat capacity cV (per unit volume), so that „e = cV „T . On the other hand, experimen- tally it is known that the energy (heat) current density is proportional to the temperature gradient. For isotropic materials this relation has the simplest form q = −k ∇T, where k is thermal conductivity depending on the material. Putting all together and neglecting the temperature dependence of k (that is justified if changes ot the temperature are small) one obtains the thermal diffusion (or heat conduction) equation ∂T ∂t  κ∆T, κ ≡ k cV , where k is thermal diffusivity having the unit m2 ës. In the stationary case this equation becomes a Laplace equation. This equation is valid in the absence of motion, otherwise one has to add the convective term v ÿ“T on the left, and then the equation for the temperature becomes coupled to the equation for the fluid dynamics above. If the thermal energy is injected into the system (e.g., as a result of burning), one has to add a corresponding source term on the right side of the equation . Physics of diffusion is similar, and the concentration of particles c satisfies the equation ∂c ∂t  κ∆c, where k is diffusivity. If particles can be carried by the media motion, the convective term v ÿ“c, similarly to the heat conduction equation. The corresponding generalization is the Smoluchowski equation. ü Heat flow between two walls As an illistration, consider the region between the y-z planes (walls) at x = 0 and x = a > 0. The temperature at the left boundary is T0 while the temperature at the right boundary is T1. If one is interested in the stationary solution of the thermal diffusion equation, one can drop the time derivative, search for the solution in the form T = T@xD and obtain the ordinary differential equation d2 T dx2  0. The solution of this equation is a linear function of x T = C1 x + C2 = T0 + HT1 − T0L x a , where the constants have been eliminated with the help of the boundary conditions. This solution corresponds to the heat flux density q = −k dT dx = −k T1 − T0 a flowing from the hot wall to the cold wall. With these boundary conditions, the stationary solution found above will be achieved from any initial condition after a sufficiently long relaxation time. One could reformulate this problem setting the heat flux density q entering the system, say, from the left and one of the temperatures, say, on the right. In this case the boundary conditions are Mathematical_physics-13-Partial_differential_equations.nb 5 dT dx @0D  − q k , T@aD = T1. The resulting solution has the form T = T1 + q k Ha − xL. ü Heat flow from a heated sphere Another example of a stationary heat flow is a sphere of radius R at temperature T0 immersed in a heat-conductive media having the temperature T¶ at infinity. Here the symmetry allows to bypass the complete solution of the heat conduction equation by using the energy conservation argument saying that in the stationary case the total outwards heat flux through any spherical surface of radius r is the same Q = 4 πr2 q = const Thus one finds q = C r2 , similarly to the electric field of a point charge found from the Gauss theorem. Now the temperature can be found integrating the equation dT dr = − q k = C kr2 , that yields T@rD = C1 + C2 r . Finding the constants from the boundary conditions one finally obtains T@rD = T∞ + HT0 − T∞L R r , similarly to the electric potential of a point charge. This solution corresponds to the total heat flux Q = −4 πr2 k dT dr = 4 πkR HT0 − T∞L. One can see that the smaller is the body, the smaller energy input is needed to support it at a temperature exceeding the temperature of the environment (a tip for building houses, to reduce heating bills). ü Maxwell equations Electromagnetic field satisfies the system of Maxwell equations. These equations become more complicated if the media strongly interacts with the electromagnetic field by being segnetoelectric (electrically polarized), ferromagnetic, or supercon- ducting. For simplicity we will consider electromagnetic field in vacuum. There are four Maxwell equations that describe generation of two different parts of the electric field and of the magnetic field by electric charges, electric currents, and temporal change of the electromagnetic field. In general, any vector field F@rD can be represented as the sum of the potential part and solenoidal part, F@rD = −∇φ + ∇ A, 6 Mathematical_physics-13-Partial_differential_equations.nb where f and A are scalar and vector potentials. The potential part has zero curl, — ä—f = 0, while the solenoidal part has zero divergence, — ÿ— äA = 0. The first Maxwell equation describes creation of the potential part of the electric field E by electric charges, ∇ ⋅ E  ρ ε0 , where r is the electric charge density and e0 is related to units. The solenoidal part of E is created by the magnetic field B changing in time, (6)∇ E  − ∂B ∂t . The next Maxwell equation states that there is no potential part in B, ∇ ⋅ B  0, because there are no magnetic charges. The last Maxwell equation describes creation of the solenoidal part of B by the electric current and electric field E changing with time, ∇ B  µ0 j + µ0 ε0 ∂E ∂t , where j is the electric current density and m0 is related to units. Representing the magnetic field in the form B = -“y + “ äA, from the third Maxwell equation one obtains Dy = 0. In the absence of the source in this equation, its solution is zero, y = 0, as said above. So that one finally has (7)B = ∇ A. Substituting this into the last Maxwell equation, one obtains the equation for the vector potential A (8)∇ ∇ A = ∇H∇ ⋅ AL − ∆A  µ0 j + µ0 ε0 ∂E ∂t . As there is a freedom in the definition of A, one can additionally require “ ÿA = 0 (the so-called calibration). After that only the vector Laplacian remains on the left side of this equation. Let us now consider equations for the electric field. The second equation becomes ∇ E  − ∂ ∂t ∇A and can be integrated (by removing “ä) to give (9)E  −∇φ − ∂A ∂t . Here the first term arises as an integration "constant" and is the potential part of the electric field. Now substituting this into the first Maxwell equation and using “·A = 0, one obtains the Poisson equation for the electric potential (10)∆φ  − ρ ε0 . Substituting Eq. (9) into Eq. (8) to eliminate E, one obtains the closed equation for the vector potential (11) 1 c2 ∂2 A ∂t2 − ∆A  µ0 j − 1 c2 ∂∇φ ∂t , c = 1 µ0 ε0 . Mathematical_physics-13-Partial_differential_equations.nb 7 Plot3D@Txt@x, tD, 8x, 0, L<, 8t, 0, tMax<, PlotRange → All, AxesLabel → AutomaticD The greater the thermal diffusivity, the faster the temperature approaches its equilibrium value. The fact that the temperature relaxes to a time-independent stationary state may be used to solve the Laplace equation and similar non-evolutionary equations by the method of relaxation. ü 1d: Temperature driven by a time-dependent boundary condition On the left side the temperature is oscillating in time and on the right side the system is thermally insulated κ = 0.1; L = 1; tMax = 20; ω = 1; Equation = ∂tT@x, tD − κ ∂x,xT@x, tD  0 BConds = 8T@0, tD  Sin@ω tD, H∂xT@x, tD ê. x → LL  0< ICond = T@x, 0D  0; Sol = NDSolve@Join@8Equation<, BConds, 8ICond<D, T, 8x, 0, L<, 8t, 0, tMax<D Txt@x_, t_D := T@x, tD ê. Sol@@1DD TH0,1L@x, tD − 0.1 TH2,0L@x, tD  0 9T@0, tD  Sin@tD, TH1,0L@1, tD  0= 88T → InterpolatingFunction@880., 1.<, 80., 20.<<, <>D<< 10 Mathematical_physics-13-Partial_differential_equations.nb Plot3D@Txt@x, tD, 8x, 0, L<, 8t, 0, tMax<, PlotRange → All, AxesLabel → Automatic, PlotPoints → 30D The greater the thermal diffusivity, the further to the right goes the disturbance. ü 2d: Solution of the Laplace equation by the relaxation method The problem is to solve the Laplace equation ∆T  0 in the rectangular 2d region 0 § x § Lx, 0 § y § Ly with the boundary conditions T@0, yD = T0, T@Lx, yD = T1, and T@x, 0D = TAx, LyE = 0. The function T can be temperature or concentration, or electric potential, or anything else that is described by the Laplace equation. The idea of the relaxation method is to solve the equation ∂T ∂t  κ∆T with some initial condition. After some time the system will come to the time independent state in which ∑t T = 0 and thus DT = 0. Setting the initial condition, one has to keep in mind that it has to be compatible with the boundary conditions, otherwise Mathematica might be unable to solve the problem. One of such initial conditions is T@x, y, 0D  T0 HLx − xL y ILy − yM x + HLx − xL y ILy − yM + T1 x y ILy − yM HLx − xL + x y ILy − yM The parameters k and tMax can be determined experimentally. In fact, one can set k = 1 and choose an appropriate tMax because these two parameters in fact enter as the product ktMax. The number of spatial discretization grid points can be controlled by the MaxSteps option. As the initial condition is singular, Mathematica uses by default 100 grid points. However, it complains anyway: Because of the singularity, any number of grid points is formally insufficient (although practically 100 grid points is OK). Apart of this, the initial condition contains division by zero and works only with an added small number in the denominator. After that initial and boudary conditions contradict each other. However, this does not really lead to a problem because the inconsistency is very small. The resulting solution is good. Mathematical_physics-13-Partial_differential_equations.nb 11 κ = 1; Lx = 1; Ly = 1; tMax = 1; T0 = 1; T1 = 0.2; Equation = ∂tT@x, y, tD − κ I∂x,xT@x, y, tD + ∂y,yT@x, y, tDM  0; BConds = 8 T@0, y, tD  T0, T@Lx, y, tD  T0, T@x, 0, tD  0, T@x, Ly, tD  0 <; ICond = T@x, y, 0D  T0 HLx − xL y HLy − yL x + HLx − xL y HLy − yL + 10−15 H∗ avoid division by 0 ∗L + T1 x y HLy − yL HLx − xL + x y HLy − yL + 10−15 H∗ avoid division by 0 ∗L ; Sol = NDSolve@Join@8Equation<, BConds, 8ICond<D, T, 8x, 0, Lx<, 8y, 0, Ly<, 8t, 0, tMax<, MaxSteps → 8110, 110, Infinity<D; Txyt@x_, y_, t_D := T@x, y, tD ê. Sol@@1DD Now check that tMax is sufficient for relaxation by plotting T at some points vs. t. Plot@8Txyt@0.5, 0.5, tD, Txyt@0.25, 0.5, tD, Txyt@0.5, 0.3, tD<, 8t, 0, tMax<, PlotRange → AllD 0.2 0.4 0.6 0.8 1.0 0.20 0.30 0.35 0.40 0.45 0.50 0.55 The relaxation is complete. One can see that the initial condition was a good approximation because the changes are not large. Now plot the initial condition tt = 0; Plot3D@Txyt@x, y, ttD, 8x, 0, Lx<, 8y, 0, Ly<, PlotRange → All, AxesLabel → AutomaticD 12 Mathematical_physics-13-Partial_differential_equations.nb Grad@f_D := 9∂xf, ∂yf= Q@x_, y_D = −Grad@Txyt@x, y, tMaxDD; H∗ Heat flux ∗L StreamPlot@8Q@x, yD@@1DD, Q@x, yD@@2DD<, 8x, 0, Lx<, 8y, 0, Ly<, VectorPoints → 10D 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 The two heat fluxes collide but the flux from the left is stronger (T0 = 1; T1 = 0.2). If we reduce the temperature of the right wall to T1 = 0.1, the flux from the left wins totally in the middle, whereas there is still a peripheric flux from the right. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 ü 1d: Deflagration (burning) Above we have considered the thermal runaway (or explosive instability) by assuming that the process is uniform in space. In reality it is not completely true because the heat exchange with the environment is not uniform and occurs at the boundary of the region with explosives. Thus the thermal instability process should begin in the center and then spread to the periph- ery. This spreading can happen in the form of propagating fronts of thermal runaway. The usual slow burning fronts that we observe igniting a sheet of paper is also described by this mechanism, it does not have to be a violent process such as explosion. In physics burning that occurs via moving fronts is called deflagration. Mathematical_physics-13-Partial_differential_equations.nb 15 Deflagration consists of two processes: (i) Decay of a flammable substance (fuel, explosive) with a rate strongly increasing with temperature and (ii) Spreading of the released heat from hot (burned) to cold (yet unburned) areas via heat conduction. The mathematical model of deflagration uses two equations: (i) Rate (or relaxational) equation for the mass of fuel m and (ii) heat conduction equation with a source m t  −Γm, Γ@TD = Γ0 ExpB− ∆U kB T F ∂T ∂t  κ∆T − Em C m t . Here Em is the energy released by burning of one unit mass of the fuel and C is the heat capacity. The expression for the burning rate of the fuel G@TD is taken to be activational (Arrhenius) with DU being the activation energy (the potential barrier to overcome for the decay). The first deflagration equation is ODE while the second one is a PDE. We will consider 1 d problem with the fuel localized in the region 0 § x § L and we can set the temperature at the ends, as well as the initial temperature, to the environmental temperature T0. Clearly T0 ` DU has to be satisfied, otherwise burning occurs immediately. Then with time the temperature in the depth of the fuel region will increase because of the field decay. If the cooling via heat conduction is sufficient, the system will reach a stationary state with the maximal temperature in the center, x = Lê2, that only slightly exceeds T0. If cooling is insufficient (k too small or L too large) the temperature will continue to increase and thermal runaway occurs. In this case burning occurs throughout the whole region practically at the same time. This is qualitatively similar to the 0 d Semenov model considered above but more involved because of the PDE. The expression for the threshold of thermal runaway in this case was obtained by Frank-Kamenetskii. If the temperature is increased or heat is injected locally, burning starts at this point and then deflagration front propagates throughout the sample. This is the most interesting case that will be modeled below. To initiate the deflagration front, we will raise the temperature at the left end of the sample, x = 0, during a short time. 16 Mathematical_physics-13-Partial_differential_equations.nb In[52]:= L = 1000; H∗ Length of the region ∗L tMax = 100; H∗ Time of calculation ∗L kB = 1; Em = 3000; H∗ Energy released per unit mass of burned fuel ∗L c = 1; H∗ Heat capacity ∗L m0 = 1; H∗ Initial mass of fuel ∗L κ = 200; H∗ Thermal diffusivity ∗L H∗ Relaxation rate ∗L Γ0 = 20; ∆U = 4000; Γ@T_D := Γ0 −∆UêT; T0 = 300; H∗ Temperature of the environment ∗L H∗ Ignition by raising temperature on the left end ∗L TIgnition = 600; tIgnition = 1; TLeft@t_D = T0 + HTIgnition − T0L TanhB t tIgnition F; H∗ Equations and solving ∗L EqT = ∂tT@x, tD + Em c ∂tm@x, tD  κ ∂x,xT@x, tD; Eqm = ∂tm@x, tD  −Γ@T@x, tDD m@x, tD; IniConds = 8T@x, 0D == T0, m@x, 0D  1<; BConds = 8T@0, tD == TLeft@tD, T@L, tD == T0<; solution = NDSolve@Join@8EqT, Eqm<, IniConds, BCondsD, 8T, m<, 8x, 0, L<, 8t, 0, tMax<, Method → 8"MethodOfLines", "SpatialDiscretization" → 8"TensorProductGrid", "MinPoints" → 200<<D mxt@x_, t_D := m@x, tD ê. First@solutionD Txt@x_, t_D := T@x, tD ê. First@solutionD Out[68]= 88T → InterpolatingFunction@880., 1000.<, 80., 100.<<, <>D, m → InterpolatingFunction@880., 1000.<, 80., 100.<<, <>D<< In[71]:= Plot3D@mxt@x, tD, 8t, 0, 0.5 tMax<, 8x, 0, L<, PlotRange → 80, 1<, PlotPoints → 100, AxesLabel → AutomaticD Out[71]= One can see a deflagration front propagating with a constant speed from the left to right end of the region. Mathematical_physics-13-Partial_differential_equations.nb 17 L = 10; λ = 20 L; H∗ Wave length. λ1=2L is the fundamental wave length in this case ∗L c = 1; H∗ Wave speed ∗L ω = 2 π c λ ; H∗ Frequency ∗L tMax = 30 L c ; H∗ Maximal time of the calculation. Lêc is the time for the wave to cross the region ∗L Equation = c2 ∂x,xφ@x, tD − ∂t,tφ@x, tD  0; H∗ 1d wave equation for the velocity potential ∗L IConds = 8φ@x, 0D  0, H∂tφ@x, tD ê. t → 0L  0<; H∗ All zero at initial time ∗L BConds = 8 H∂xφ@x, tD ê. x → 0L == Tanh@tD Sin@ ω tD, H∗ The excitation begins gradually to avoid inconsistency of initial and boundary conditions ∗L H∂xφ@x, tD ê. x → LL  0 H∗ Closed on the right end ∗L <; Timing@Sol = NDSolve@Join@8Equation<, BConds, ICondsD, φ, 8x, 0, L<, 8t, 0, tMax<, Method → 8"MethodOfLines", "SpatialDiscretization" → 8"TensorProductGrid", "MinPoints" → 100<<D;D φxt@x_, t_D := φ@x, tD ê. Sol@@1DD vxt@x_, t_D := −∂xxφxt@xx, tD ê. xx → x H∗ Velocity of the media in the wave is negative gradient of its potential ∗L 87.188, Null< Calculation with λ = 23ê2 L - longer than the fundamental wave length, nonresonant case. The wave pattern is irregular and the wave amplitude does not increase with time, remaining of the order of the excitation amplitude at the left end. Plot3D@vxt@x, tD, 8x, 0, L<, 8t, 0, tMax<, PlotPoints → 100, AxesLabel → AutomaticD Calculation with λ = 2 L - the fundamental wave length, resonance. The wave pattern is regular with nodes at both bound- aries, and the wave amplitude unlimitedly increases with time, as always the case for resonance in linear undamped systems. 20 Mathematical_physics-13-Partial_differential_equations.nb Plot3D@vxt@x, tD, 8x, 0, L<, 8t, 0, tMax<, PlotPoints → 100, AxesLabel → AutomaticD Calculation with λ = L - the second overtone. The wave pattern is regular with nodes at both boundaries and in the middle, and the wave amplitude unlimitedly increases with time. Plot3D@vxt@x, tD, 8x, 0, L<, 8t, 0, tMax<, PlotPoints → 100, AxesLabel → AutomaticD Let us now consider excitation of standing waves in a pipe with the left end exciting waves as before and the right end open. As before, we consider the problem as 1d neglecting the influence of side walls. The boundary condition at the right end changes to f@L, t] ã 0. To see it, use Eq. (5), ∑tf ∂ dP. At the open end, the pressure with a good accuracy (although not absolutely) merges with the constant atmospheric pressure, thus dP@L, tDã 0. Thus also ∑tf@L, tD ã 0. Integrating over time, one obtains f@L, tDã const. However, the potential is defined up to an arbitrary constant, so that this constant can be set to zero. The results show that the wave amplitude unlimitedly increases with time if the wave length coincides with the fundamental wave length Mathematical_physics-13-Partial_differential_equations.nb 21 λ1 = 4 L or those of the overtones, ln = l1 n , n = 1, 3, 5, — Heven overtones absentL This is a resonance behavior. L = 10; λ = 4 L; H∗ Wave length. λ1=4L is the fundamental wave length in this case ∗L c = 1; H∗ Wave speed ∗L ω = 2 π c λ ; H∗ Frequency ∗L tMax = 30 L c ; H∗ Maximal time of the calculation. Lêc is the time for the wave to cross the region ∗L Equation = c2 ∂x,xφ@x, tD − ∂t,tφ@x, tD  0; H∗ 1d wave equation for the velocity potential ∗L IConds = 8φ@x, 0D  0, H∂tφ@x, tD ê. t → 0L  0<; H∗ All zero at initial time ∗L BConds = 8 H∂xφ@x, tD ê. x → 0L == Tanh@tD Sin@ ω tD, H∗ The excitation begins gradually to avoid inconsistency of initial and boundary conditions ∗L φ@L, tD  0 H∗ Closed on the right end ∗L <; Timing@Sol = NDSolve@Join@8Equation<, BConds, ICondsD, φ, 8x, 0, L<, 8t, 0, tMax<, Method → 8"MethodOfLines", "SpatialDiscretization" → 8"TensorProductGrid", "MinPoints" → 100<<D;D φxt@x_, t_D := φ@x, tD ê. Sol@@1DD vxt@x_, t_D := −∂xxφxt@xx, tD ê. xx → x H∗ Velocity of the media in the wave is negative gradient of its potential ∗L 82.61, Null< Calculation with λ = 4 L - the fundamental. The wave pattern is regular with a node at the left and an anti-node at the right ends of the pipe, and the wave amplitude unlimitedly increases with time. Plot3D@vxt@x, tD, 8x, 0, L<, 8t, 0, tMax<, PlotPoints → 100, AxesLabel → AutomaticD 22 Mathematical_physics-13-Partial_differential_equations.nb tt = tMax; Plot3DBφxyt@x, y, ttD, 8x, −Lx, Lx<, 8y, 0, Ly<, PlotRange → All, AxesLabel → Automatic, PlotPoints → 100, AspectRatio → Ly 2 Lx F Plotting the wave pattern together with the theoretical result for the directions of the minima Mathematical_physics-13-Partial_differential_equations.nb 25 tt = 10; ShowB DensityPlotBφxyt@x, y, ttD, 8x, −Lx, Lx<, 8y, 0, Ly<, PlotRange → All, PlotPoints → 100, AspectRatio → Ly 2 Lx F, PlotB:CotBArcSinB λ d0 FF x, −CotBArcSinB λ d0 FF x, CotBArcSinB 2 λ d0 FF x, −CotBArcSinB 2 λ d0 FF x>, 8x, −10, 10<, PlotRange → 88−Ly, Ly<, 80, Lx<<, AspectRatio → Lx 2 Ly F F Let us create an animated GIF file to put it onto a web site. 26 Mathematical_physics-13-Partial_differential_equations.nb H∗ Preparing to export into animated GIF ∗L H∗ Set PlotRange Hto avoid jitterL, Viewpoint, BoxRatios, and ImageSize ∗L tt = tMax; Plot3DBφxyt@x, y, ttD, 8x, −Lx, Lx<, 8y, 0, Ly<, PlotRange → 8−0.26, 0.26<, PlotPoints → 100, BoxRatios → :1, Ly 2 Lx , 0.3>, ViewPoint → 80, −100, 100<, ImageSize → 500F H∗ Export into animated GIF ∗L H∗ Create a list of frames, Plot3Ds at different times ∗L WaveFrame@t_D := Plot3DBφxyt@x, y, tD, 8x, −Lx, Lx<, 8y, 0, Ly<, PlotRange → 8−0.26, 0.26<, PlotPoints → 100, BoxRatios → :1, Ly 2 Lx , 0.3>, ViewPoint → 80, −100, 100<, Axes → None, ImageSize → 500F NFrames = 100; tFrame@k_D := tMax k NFrames ; WaveFrames = Table@WaveFrame@tFrame@kDD, 8k, 1, NFrames<D; Export@"D:\\_Lehman\\2010 Fall\\PHY 307\\Notes\Diffraction−one−slit.gif", WaveFramesH∗,"TransparentColor"→White∗LD ; ü 2d: Two-slit diffraction in time Two-slit diffraction results in the maxima of the wave intensity for the directions specified by the angle q that satisfy d Sin@θD = m λ, m = 0, ±1, ±2, — Here d is the distance between the slits and l is the wave length. We solve the wave equation in a rectangular region with the slips in the middle of the bottom line within the time interval between opening the slits and the wave hitting the far (top) wall. (Stationary wave equation cannot be solved by Mathematica at the moment.) Mathematical_physics-13-Partial_differential_equations.nb 27 tt = 10; ShowB DensityPlotBφxyt@x, y, ttD, 8x, −Lx, Lx<, 8y, 0, Ly<, PlotRange → All, PlotPoints → 100, AspectRatio → Ly 2 Lx F, PlotB:1010 x, CotBArcSinB λ d FF x, −CotBArcSinB λ d FF x, CotBArcSinB 2 λ d FF x, −CotBArcSinB 2 λ d FF x>, 8x, −10, 10<, PlotRange → 88−Ly, Ly<, 80, Lx<<, AspectRatio → Lx 2 Ly F F ü Schrödinger equation ü 2d: One-slit diffraction in time Since the original complex Schrödinger equation cannot be solved by NDSolve, its real and complex parts have to be separated that leads to doubling the equation and boundary and initial conditions. We consider a plane wave, Eq. (14) falling on a slit from y  0 and solve the time-dependent Schrödinger equation, Eq. (12) on the other side of the slit in the region 0 § y § Ly and -Lx § x § Lx. We assume that the region is bounded by rigid walls, so that the boundary condition is Yã 0 everywhere except the slit of width d0 in the bottom wall, where is Y is the same as in the incident plane wave, Y~Exp@-ÂwtD. It is convenient to express the frequency w through the particle's de Broglie wave length l as ω = — 2 m 2 π λ 2 , see Eqs. (15) and (16). Having l as a parameter in calculations, we can relate numerical result to the one-slit diffraction formula, Eq. (20). We gradually open the slit starting at t = 0, to avoid inconsistency of the initial and boundary conditions. As the initial condition we take Y@x, y, 0D ã 0, no particles in the box. 30 Mathematical_physics-13-Partial_differential_equations.nb Lx = 10; Ly = 10; H∗ Sizes of the rectangular spatial region ∗L — = 1; H∗ Planck's constant ∗L m = 1; H∗ Mass of the particle ∗L d0 = 3; H∗ Width of the slit ∗L λ = 1; H∗ De Broglie wave length of the particle ∗L ω = — 2 m 2 π λ 2 ; H∗ Frequency corresponding to de Broglie wave length ∗L tMax = 3; H∗ Time of the calculation ∗L H∗ Modeling the slit ∗L PowSlit = 50; FSlit@x_D = 1 1 + H2 x ê d0L2 PowSlit ; Plot@FSlit@xD, 8x, −Lx, Lx<, PlotRange → AllD tOpenSlit = 0.1; H∗ Time to fully open the slit ∗L H∗ Creating and solving the equations for the velocity potential ∗L H∗ Boundary conditions: perpendicular component of the velocity on boundaries is zero, except for the slit and is a given periodic function on the slit ∗L Equations = : — ∂tΨRe@x, y, tD + —2 2 m I∂x,xΨIm@x, y, tD + ∂y,yΨIm@x, y, tDM  0, −— ∂tΨIm@x, y, tD + —2 2 m I∂x,xΨRe@x, y, tD + ∂y,yΨRe@x, y, tDM  0 >; IConds = 8ΨRe@x, y, 0D  0, ΨIm@x, y, 0D  0<; H∗ All zero at initial time ∗L BConds = : ΨRe@x, 0, tD == FSlit@xD TanhB t tOpenSlit F Cos@ω tD, ΨIm@x, 0, tD  −FSlit@xD TanhB t tOpenSlit F Sin@ω tD, H∗ The slit opens gradually to avoid inconsistency of initial and boundary conditions ∗L ΨRe@x, Ly, tD  0, ΨIm@x, Ly, tD  0, ΨRe@Lx, y, tD  0, ΨIm@Lx, y, tD  0, ΨRe@−Lx, y, tD  0, ΨIm@−Lx, y, tD  0 >; Timing@Sol = NDSolve@Join@Equations, BConds, ICondsD, 8ΨRe, ΨIm<, 8x, −Lx, Lx<, 8y, 0, Ly<, 8t, 0, tMax<, Method → 8"MethodOfLines", "SpatialDiscretization" → 8"TensorProductGrid", "MinPoints" → 80<<D;D ΨRexyt@x_, y_, t_D := ΨRe@x, y, tD ê. Sol@@1DD ΨImxyt@x_, y_, t_D := ΨIm@x, y, tD ê. Sol@@1DD Ψ2xyt@x_, y_, t_D := ΨRexyt@x, y, tD2 + ΨImxyt@x, y, tD2 Mathematical_physics-13-Partial_differential_equations.nb 31 -10 -5 5 10 0.2 0.4 0.6 0.8 1.0 8289.36, Null< tt = tMax; Plot3DBΨ2xyt@x, y, ttD, 8x, −Lx, Lx<, 8y, 0, Ly<, PlotRange → 80, 2<, PlotPoints → 100, BoxRatios → :1, Ly 2 Lx , 0.3>, ViewPoint → 80, −100, 100<, ImageSize → 500F The results show that quantum particles entering the box via the slit undergo diffraction with the directions of diffraction minima in accordance with the elementary wave theory. Before the quantum wave hits the far wall, the plotted probability density †Y@r, tD§2 remains smooth. However, after the beam of particles hits the far wall, because of the interference of the incident and reflected wave the wave pattern becomes oscillating in space. The total probability for a particle to be in the box Ÿ †Y@r, tD§2 „r steadily increases as more and more particles are entering the box via the slit. Let us create an animated GIF file to put it onto a web site. 32 Mathematical_physics-13-Partial_differential_equations.nb