Partial differential equations

This chapter is an introduction to PDE with physical examples that allow straightforward numerical solution with Mathematica. 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,…]==0

(Δ is the Laplace operator), Poisson equation (Laplace equation with a source)

Δϕ[x,y,…]==f[x,y,…],

wave equation

(∂_t^2)ϕ[t,x,y,…]-(c^2)Δϕ[t,x,y,…]==0,

heat conduction / diffusion equation

(∂_t)ϕ[t,x,y,…]-κΔϕ[t,x,y,…]==0,

Schrödinger equation

i(∂_t)ϕ[t,x,y,…]+(aΔ+bf[x,y,…])ϕ[t,x,y,…]==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,y))f[x,y]==0

is

f[x,y]=F[x]+G[y],

where F[x] and G[y] are arbitrary functions. The solution of the first-order PDE

(∂_t)f[t,x]-v(∂_x)f[t,x]==0

is

f[t,x]=g[x-vt]

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. Mathematica 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, electromagnetic fields, and quantum-mechanical probabilities. Below we review these equations.

Continuity equation

Any conserved quantity described by the density ρ and the corresponding current j satisfies the continuity equation

((∂ρ)/(∂t))+divj==0.

This is the continuity equation in the differential form. If dV is an infinitesimal volume, (∂_t)ρdV is the rate of change of the amount (corresponding to ρ) within this volume, whereas divj dV 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))∫ρdV==- ∫divjdV.

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))∫ρdV==- ∫j·dS.

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 (), there ρ is the mass density and the mass current j density is given by

j=ρv.

Liquids, unlike gases, are practically uncompressible, ρ=const, thus the continuity equation simplifies to

divv==0.

Dynamics of a non-viscous fluid is described by the Euler equation

((∂v)/(∂t))+(v∇)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 (v∇)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

((∂v)/(∂t))+(v∇)v==-((∇P)/ρ)+g+(η/ρ)Δv,

where η is the viscosity coefficient.

Euler equation has a limited applicability because its solutions (for instance, for an object in a flow) have discontinuities of the velocity in the stream. In these regions the viscosity becomes important and at high enough speeds the flow becomes turbulent. Still, Euler equation generates Bernoilli equation that explains the lifting force on airplanes as well as the mechanism of sound creation by voice bands.

Navier-Stokes equation is more realistic and it can be applied, for instance, to the problem of a drag force acting on a object moving in a liquid. For small boby speeds u the quadratic terms can be dropped, equation becomes linear and can be solved analytically for simple shapes such as sphere. The drag forces on a solid sphere and air bubble of radius R have the form

F=(Piecewise[{{6πηRu,, solid sphere}, {4πηRu,, air bubble.}}])

The difference in the coefficient is due to different bounary conditions: The liquid sticks to the solid sphere but glides along the surface of the air bubble that reduces the resistance.

Waves

Waves exist in compressible media such as fluids and solids. The main prerequisite for the waves' existence is a restoring force. In fluids, there is only a restoring force with respect to change of the volume, thus the waves are waves of compression / expansion. In these waves the media displacement is collinear with the direction of the wave propagation, so that these waves are called longitudinal waves. In solids, there is additionally a restoring force with respect to change of the shape of a media element. The corresponding additional waves are transverse waves. Below we will consider small-amplitude longitudinal waves in fluids.

In small-amplitude waves pressure and density only slightly deviate from their equilibrium values, so that one can define

P=(Subscript[P, 0])+δP,        ρ=(Subscript[ρ, 0])+δρ

and linearize equations in small δP and δρ, neglecting quadratic velocity terms because the velocity is small. The continuity equation becomes

(∂/(∂t))δρ+(Subscript[ρ, 0])divv==0,

whereas in the absence of viscosity from the Euler equation one obtains

((∂v)/(∂t))==-((∇δP)/(Subscript[ρ, 0])).

The small changes of pressure and density are related by the adiabatic compressibility,

δP=(Subscript[(((∂P)/(∂ρ))), S])δρ.

One can obtain a single equation describing waves by differentiating the continuity equation over time and then using the Euler equation. This yields the wave equation

(∂_t^2)δρ-(c^2)Δδρ,           c≡(Sqrt[(Subscript[(((∂P)/(∂ρ))), S])]),

where c is the speed of sound. Note that the speed of sound (that can be large) has no relation to the velocity of the media (that is small). Since δP~δρ, it satisfies the same equation,

(∂_t^2)δP-(c^2)ΔδP.

Let us obtain the equation for the velocity. Taking rotor of the linearized Euler equation, one obtains

(∂/(∂t))∇×v==0,

thus velocity is a potential field, ∇×v=0, and it can be searched in the form

v=-∇ϕ.

Inserting it into the linearized Euler equation () and integrating it by removing the gradient one obtains

((∂ϕ)/(∂t))==((δP)/(Subscript[ρ, 0]))=((c^2)/(Subscript[ρ, 0]))δρ.

Differentiating this over time and using the continuity equation, one obtains the same wave equation for the velocity potential,

(∂_t^2)ϕ-(c^2)Δϕ.

Let us now seek for the solution of the wave equation in the form of a plane wave

ϕ=ACos[ωt-k·r+(Subscript[ϕ, 0])],

where A is the amplitude, ω is the frequency and k is the wave vector that defines the direction of the wave propagation. After substitution, differentiation, and cancellation one obtains the dispersion relation

(ω^2)-(c^2)(k^2)==0.

The period T and wave length λ are related to ω and k by

ω=((2π)/T),     k=((2π)/λ).

One can see that ω can be different while k adjusts to ω as ck=ω.

Let us now find the media velocity in a plane wave:

v=-∇ϕ=-A k Sin[ωt-k·r+(Subscript[ϕ, 0])].

One can see that v is collinear to k, this is why such waves are called longitudinal waves.

Since the wave equation is linear, a superposition of plane waves moving in different directions with different frequencies and wave lengths is also its solution. This illustrates the fact that there is much freedom in the solution of PDEs. The selection of the actual solution is done by the initial and boundary conditions.

In the important case of monochromatic waves whose time dependence is describe by a single frequency,

ϕ[r,t]=ϕ[r]Sin[ωt],

or similar with Sin[ωt] or Exp[iωt], the wave equation simplifies to the time-independent stationary wave equation

Δϕ+(k^2)ϕ==0,                k=(ω/c).

In one dimension it becomes an ODE.

Thermal conductivity and diffusion

Temperature in an isolated body tends to equilibrate with time because the energy flows from warmer regions to colder regions. Since the energy is conserved, its flow obeys the continuity equation

((∂e)/(∂t))+div q==0,

where e is the density of the internal energy and q≡(Subscript[j, e]) is its current. If the body is incompressible, the change of the internal energy according to thermodynamics is due to the heat received, de=δQ. On the other hand, the heat received is related to the temperature change via heat capacity (Subscript[c, V]) (per unit volume), so that de=(Subscript[c, V])dT. On the other hand, experimentally 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/(Subscript[c, V])),

where κ is thermal diffusivity having the unit (m^2)/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 κ 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 (Subscript[T, 0]) while the temperature at the right boundary is (Subscript[T, 1]). 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[x] and obtain the ordinary differential equation

(((d^2)T)/(dx^2))==0.

The solution of this equation is a linear function of x

T=(Subscript[C, 1])x+(Subscript[C, 2])=(Subscript[T, 0])+((Subscript[T, 1])-(Subscript[T, 0]))(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(((Subscript[T, 1])-(Subscript[T, 0]))/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

(dT/dx)[0]==-(q/k),         T[a]=(Subscript[T, 1]).

The resulting solution has the form

T=(Subscript[T, 1])+(q/k)(a-x).

Heat flow from a heated sphere

Another example of a stationary heat flow is a sphere of radius R at temperature (Subscript[T, 0]) immersed in a heat-conductive media having the temperature (Subscript[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((πr)^2)q=const

Thus one finds

q=(C/(r^2)),

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/(kr^2)),

that yields

T[r]=(Subscript[C, 1])+((Subscript[C, 2])/r).

Finding the constants from the boundary conditions one finally obtains

T[r]=(Subscript[T, ∞])+((Subscript[T, 0])-(Subscript[T, ∞]))(R/r),

similarly to the electric potential of a point charge. This solution corresponds to the total heat flux

Q=-4((πr)^2)k(dT/dr)=4πkR((Subscript[T, 0])-(Subscript[T, ∞])).

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 superconducting. 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[r] can be represented as the sum of the potential part and solenoidal part,

F[r]=-∇ϕ+∇×A,

where ϕ and A are scalar and vector potentials. The potential part has zero curl, ∇×∇ϕ=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==(ρ/(Subscript[ϵ, 0])),

where ρ is the electric charge density and (Subscript[ϵ, 0]) is related to units. The solenoidal part of E is created by the magnetic field B changing in time,

∇×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==(Subscript[μ, 0])j+(Subscript[μ, 0])(Subscript[ϵ, 0])((∂E)/(∂t)),

where j is the electric current density and (Subscript[μ, 0]) is related to units.

Representing the magnetic field in the form B=-∇ψ+∇×A, from the third Maxwell equation one obtains Δψ=0. In the absence of the source in this equation, its solution is zero, ψ=0, as said above. So that one finally has

B=∇×A.

Substituting this into the last Maxwell equation, one obtains the equation for the vector potential A

∇×∇×A=∇(∇·A)-ΔA==(Subscript[μ, 0])j+(Subscript[μ, 0])(Subscript[ϵ, 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

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

Δϕ==-(ρ/(Subscript[ϵ, 0])).

Substituting Eq. () into Eq. () to eliminate E, one obtains the closed equation for the vector potential

(1/(c^2))(((∂^2)A)/(∂(t^2)))-ΔA==(Subscript[μ, 0])j-(1/(c^2))((∂∇ϕ)/(∂t)),          c=(1/(Sqrt[(Subscript[μ, 0])(Subscript[ϵ, 0])])).

In the static case time derivatives vanishes and this equation describes generation of the vector potential A and thus of the magnetic field B by the electric current. In general, this is a wave equation for electromagnetic waves with c being the speed of light. There are two sources on the right side of this equation that excite electromagnetic waves. One source is electric current and another is time-dependent potential part of the electirc field. The latter is created by moving electric charges according to Eq. ().

As above, one can seek for the solution of this wave equation in the form of a plane wave. It turns out that in this plane wave E and B are perpendicular to each other and to the wave vector k. This means that electromagnetic waves are transverse waves (Excercise).

Since the electric current is due to the motion of electric changes in space, one can write the continuity equation () with the electric current density j=ρv.

Schrödinger equation

In quantum mechanics, dynamics of an isolated system is described by the Schrödinger equation for the complex wave function Ψ[r,t] that has the form

iℏ((∂Ψ)/(∂t))==(OverHat[H])Ψ,

where ℏ is the Plank constant and (OverHat[H]) is the so-called Hamilton operator, the energy of the system expressed via positions and momenta, whereas the momenta are replaced by first-order spatial differential operators. For instance, for a particle with a potential energy U[r] the Hamiltonian function (the energy) is given by

H=((p^2)/(2m))+U[r].

Replacing p by the quantum-mechanical operator

(OverHat[p])=-iℏ∇,

one obtains the Hamilton operator

(OverHat[H])=-(((ℏ^2)Δ)/(2m))+U[r]

that acts on the wave function in the Schrödinger equation. This equation resembles the heat conduction and diffusion equations, only it has an imaginary factor at the time derivative. In spite of this similarity, there is little in common between the two equations.Whereas the heat conduction equation is relaxational, Schrödinger equation describes a conservative dynamics.

The interpretation of the wave function is that ((|Ψ[r, t]|)^2) ≡ Ψ[r, t] (SuperStar[Ψ[r, t]]) is the probability density for a particle to be at the position r. Practically this is not of the primary importance, however. What mostly matters in quantum mechanics is stationary states that depend of time monochromatically,

Ψ[r,t]=ψ[r]Exp[((iE)/ℏ)t],

E being the energy of the particle, and satisfy the stationary Schrödinger equation

(OverHat[H])ψ==Eψ.

These states corresponds to quantum energy levels (next chapter). In fact, Schrödinger equation has more in common with the wave equation and it also has a solution in the form of a plane wave if U[r]=U=const

Ψ[r,t]=A Exp[-i(ωt-k·r)].

Substituting this into Schrödinger equation one obtains the relation

ℏω=(((ℏ^2)(k^2))/(2m))+U

between the frequency ω and the wave vector k. Frequency ω is related to the energy E and the wave vector to the de Broglie wave length as

ℏω=E,     k=2π/λ.

Solving Eq. () for k one obtains

k=(Sqrt[2m(E-U)])/ℏ.

E>U is the regular case in which the kinetic energy E-U is positive and k is real. The case E<U corresponds to the negative kinetic energy and describes quantum motion in the regions where classical motion is impossible (e.g., under a potential barrier). In this case k is imaginary,

k=±iκ,          κ=(Sqrt[2m(U-E)])/ℏ,

and the plane wave is an exponential imaginary in time but real in space:

Ψ[r,t]=A Exp[-iωt+κ·r].

There are a lot of different forms of the Schrödinger equation. For instance, in the case of rotational motion the angular momentum becomes a differential operator. For a system of interacting particles the wave function depends on all positions and there are Laplacians for each particle. This multidimensional PDE cannot be solved directly without drastic approximations, analytically and numerically.

Quantum mechanics of relativistic particles is described by Dirac equation that contains Schrödinger equation as its non-relativistic limit.

Numerical solution of PDE

Heat conduction and Laplace equations

1d: Relaxation of temperature from an initial state

It is important that the initial condition does not contradict the boundary conditions, otherwise some conditions will be discarded.

κ=0.03;    L=1;    tMax=1;
Equation=(∂_t)T[x,t]-κ(∂_(x,x))T[x,t]==0
BConds={T[0,t]==0,T[L,t]==0};
ICond=T[x,0]==Sin[((5π x)/L)];

Sol=NDSolve[Join[{Equation},BConds,{ICond}],T,{x,0,L},{t,0,tMax}]
Txt[x_,t_]:=T[x,t]/.Sol[[1]]

Plot3D[Txt[x,t],{x,0,L},{t,0,tMax},PlotRange->All,AxesLabel->Automatic]

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=(∂_t)T[x,t]-κ(∂_(x,x))T[x,t]==0
BConds={T[0,t]==Sin[ω t],((∂_x)T[x,t]/.x->L)==0}
ICond=T[x,0]==0;

Sol=NDSolve[Join[{Equation},BConds,{ICond}],T,{x,0,L},{t,0,tMax}]
Txt[x_,t_]:=T[x,t]/.Sol[[1]]

Plot3D[Txt[x,t],{x,0,L},{t,0,tMax},PlotRange->All,AxesLabel->Automatic,PlotPoints->30]

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 ≤ (Subscript[L, x]), 0≤y≤(Subscript[L, y]) with the boundary conditions T[0,y]=(Subscript[T, 0]), T[(Subscript[L, x]),y]=(Subscript[T, 1]), and T[x,0]=T[x,(Subscript[L, y])]=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 ΔT = 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,0]==(Subscript[T, 0])((((Subscript[L, x])-x)y((Subscript[L, y])-y))/(x+((Subscript[L, x])-x)y((Subscript[L, y])-y)))+(Subscript[T, 1])((x y((Subscript[L, y])-y))/(((Subscript[L, x])-x)+x y((Subscript[L, y])-y)))

The parameters κ and tMax can be determined experimentally. In fact, one can set κ=1 and choose an appropriate tMax because these two parameters in fact enter as the product κtMax.

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.

κ=1;    Lx=1; Ly=1;   tMax=1;   (Subscript[T, 0])=1; (Subscript[T, 1])=0.2;
Equation=(∂_t)T[x,y,t]-κ((∂_(x,x))T[x,y,t]+(∂_(y,y))T[x,y,t])==0;
BConds={   T[0,y,t]==(Subscript[T, 0]),
                    T[Lx,y,t]==(Subscript[T, 0]),
                    T[x,0,t]==0,
                    T[x,Ly,t]==0       };
ICond=T[x,y,0]==(Subscript[T, 0])(((Lx-x)y(Ly-y))/(x+(Lx-x)y(Ly-y)+(10^(-15)) (* avoid division by 0 *)))+(Subscript[T, 1])((x y(Ly-y))/((Lx-x)+x y(Ly-y)+(10^(-15)) (* avoid division by 0 *)));

Sol=NDSolve[Join[{Equation},BConds,{ICond}],T,{x,0,Lx},{y,0,Ly},{t,0,tMax},MaxSteps->{110,110,Infinity}];
Txyt[x_,y_,t_]:=T[x,y,t]/.Sol[[1]]

Now check that tMax is sufficient for relaxation by plotting T at some points vs. t.

Plot[{Txyt[0.5,0.5,t],Txyt[0.25,0.5,t],Txyt[0.5,0.3,t]},{t,0,tMax},PlotRange->All]

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,tt],{x,0,Lx},{y,0,Ly},PlotRange->All,AxesLabel->Automatic]

Now plot the final solution

tt=tMax;
Plot3D[Txyt[x,y,tt],{x,0,Lx},{y,0,Ly},PlotRange->All,AxesLabel->Automatic]

One can see that the final solution does not differ much from the initial condition.

There is a more elegant version of the method of relaxation that avoids the thorny problem of constructing an initial condition that satisfies the boundary conditions. The idea is to start with trivial boundary and initial conditions and then transform the boundary conditions into true ones, as is implemented in the code below.

There is a problem, however, that Mathematica does not detect from the initial and bondary conditions that there is a singularity in the problem and chooses a grid with a too small number of points. To correct it, one has to find right options in the Mathematica help. As a result, Mathematica complains because of the singularity but the solution is good.

κ=1;    Lx=1; Ly=1;   tMax=1;   (Subscript[T, 0])=1;  (Subscript[T, 1])=0.2;
Equation=(∂_t)T[x,y,t]-κ((∂_(x,x))T[x,y,t]+(∂_(y,y))T[x,y,t])==0;
BConds={    T[0,y,t]==(Subscript[T, 0])Tanh[10t],
                    T[Lx,y,t]==(Subscript[T, 1])Tanh[10t],
                    T[x,0,t]==0,
                    T[x,Ly,t]==0       };
ICond=T[x,y,0]==0;

Sol=NDSolve[Join[{Equation},BConds,{ICond}],T,{x,0,Lx},{y,0,Ly},{t,0,tMax}, Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid","MinPoints"->100}}];
Txyt[x_,y_,t_]:=T[x,y,t]/.Sol[[1]]

There are less complaints than in the main version of the relaxation method. Now check that tMax is sufficient for relaxation by plotting T at some points vs. t.

Plot[{Txyt[0.5,0.5,t],Txyt[0.25,0.5,t],Txyt[0.5,0.3,t]},{t,0,tMax},PlotRange->All]

Relaxation is complete.

tt=tMax;
Plot3D[Txyt[x,y,tt],{x,0,Lx},{y,0,Ly},PlotRange->All,AxesLabel->Automatic]

Let us now plot the heat flux (or the electric field in the case of the electric potential)

Grad[f_]:={(∂_x)f,(∂_y)f}
Q[x_,y_]=-Grad[Txyt[x,y,tMax]];  (* Heat flux *)
StreamPlot[{Q[x,y][[1]],Q[x,y][[2]]},{x,0,Lx},{y,0,Ly},VectorPoints->10]

The two heat fluxes collide but the flux from the left is stronger ((Subscript[T, 0])=1;  (Subscript[T, 1])=0.2). If we reduce the temperature of the right wall to (Subscript[T, 1])=0.1, the flux from the left wins totally in the middle, whereas there is still a peripheric flux from the right.

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 periphery. 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.

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

((dm)/(dt))==-Γm,                            Γ[T]=(Subscript[Γ, 0])Exp[-((ΔU)/((Subscript[k, B])T))]
((∂T)/(∂t))==κΔT-((Subscript[E, m])/C)((dm)/(dt)) .

Here (Subscript[E, m]) 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 Γ[T] is taken to be activational (Arrhenius) with ΔU 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 1d 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 (Subscript[T, 0]). Clearly (Subscript[T, 0])≪ΔU 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 (Subscript[T, 0]). If cooling is insufficient (κ 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 0d 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.

L=1000;  (* Length of the region *)
tMax=100;  (* Time of calculation *)
kB=1;
Em=3000;   (* Energy released per unit mass of burned fuel *)
c=1;       (* Heat capacity *)  
m0=1;     (* Initial mass of fuel *)      
κ=200;  (* Thermal diffusivity *)

(* Relaxation rate *) 
Γ0=20;      ΔU=4000;
Γ[T_]:=Γ0 (E^(-ΔU/T));

T0=300;   (* Temperature of the environment *)

(* Ignition by raising temperature on the left end *)
TIgnition=600;     tIgnition=1;
TLeft[t_]=T0+(TIgnition-T0) Tanh[(t/tIgnition)];

(* Equations and solving *)
EqT=(∂_t)T[x,t]+(Em/c)(∂_t)m[x,t]==κ(∂_(x,x))T[x,t];
Eqm=(∂_t)m[x,t]==-Γ[T[x,t]]m[x,t];
IniConds={T[x,0]==T0,m[x,0]==1};
BConds={T[0,t]==TLeft[t],T[L,t]==T0};
solution=NDSolve[Join[{EqT,Eqm},IniConds,BConds],{T,m},{x,0,L},{t,0,tMax},Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid","MinPoints"->200}}]
mxt[x_,t_]:=m[x,t]/.First[solution]
Txt[x_,t_]:=T[x,t]/.First[solution]

Plot3D[mxt[x,t],{t,0,0.5tMax},{x,0,L},PlotRange->{0,1},PlotPoints->100,AxesLabel->Automatic]

One can see a deflagration front propagating with a constant speed from the left to right end of the region.

Plot3D[Txt[x,t],{t,0,tMax},{x,0,L},PlotRange->All,PlotPoints->100,AxesLabel->Automatic]

Temperature sharply increases as a result of burning but then begins to decrease because of the energy loss to the environment, especially at the ends. Let us plot spatial profiles of the amount of fuel and temperature at a particular time.

tt=30;
Plot[mxt[x,tt],{x,0,L}]

tt=30;
Plot[Txt[x,tt],{x,0,L}]

Wave equation

Stationary wave equation for monochromatic waves Δϕ+(k^2)ϕ==0 cannot be solved directly by Mathematica's NDSolve and, unlike the Laplace equation, it cannot be solved by the handy relaxation method, also using NDSolve. Solution of the stationary wave equation requires special methods such as spatial discretization by the user that reduces it to a system of linear algebraic equations that can be then solved numerically by Mathematica. Below we will consider only the full time-dependent wave equation with initial conditions (an evolution equation) that can be solved by NDSolve.

1d: Excitation of standing waves in a closed region

Consider a closed 1d region between two walls, 0≤x≤L. The left wall makes a small-amplitude sinusoidal motion with a frequency ω that excites waves in the system. These waves reflect from the right wall that leads to formation of standing waves. The displacement of the left wall is assumed to be much smaller than L and the wave length λ=2πc/ω, where c is the speed of the wave. Thus this displacement can be ignored and we set a boundary condition at x=0. On the other hand, motion of the left wall generates motion of the media, so that the boundary condition at x=0 reads v[0,t]=(Subscript[v, wall])[t]. To avoid the problem of inconsistency of initial and boundary conditions, we switch on oscillations of the left wall gradually from zero. We solve the wave equation for the media velocity potential ϕ and them find the media velocity as v[x,t] =-(∂_x)ϕ[x,t]. This model describes standing waves in a pipe with both ends closed, whereas the influence of the side walls is neglected.

The results show that the wave amplitude unlimitedly increases with time if the wave length coincides with the fundamental wave length

(Subscript[λ, 1])=2L

or those of the overtones,

(Subscript[λ, n])=((Subscript[λ, 1])/n),     n=1,2,3,⃛

This is a resonance behavior.

L=10;
λ=(2^0)L; (* Wave length.  (Subscript[λ, 1])=2L is the fundamental wave length in this case *)
c=1;     (* Wave speed *)
ω=((2π c)/λ);   (* Frequency *)
tMax=30(L/c);   (* Maximal time of the calculation. L/c is the time for the wave to cross the region *)

Equation=(c^2)(∂_(x,x))ϕ[x,t]-(∂_(t,t))ϕ[x,t]==0;    (* 1d wave equation for the velocity potential *)
IConds={ϕ[x,0]==0,((∂_t)ϕ[x,t]/.t->0)==0};    (* All zero at initial time *)
BConds={  ((∂_x)ϕ[x,t]/.x->0)==Tanh[t]Sin[ ω t], (* The excitation begins gradually to avoid inconsistency of initial and boundary conditions *)
                  ((∂_x)ϕ[x,t]/.x->L)==0      (* Closed on the right end *)
      };

Timing[Sol=NDSolve[Join[{Equation},BConds,IConds],ϕ,{x,0,L},{t,0,tMax}, Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid","MinPoints"->100}}];]
ϕxt[x_,t_]:=ϕ[x,t]/.Sol[[1]]
vxt[x_,t_]:=-(∂_xx)ϕxt[xx,t]/.xx->x   (* Velocity of the media in the wave is negative gradient of its potential *)

Calculation with λ=(2^(3/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,t],{x,0,L},{t,0,tMax},PlotPoints->100,AxesLabel->Automatic]

Calculation with λ=2L - the fundamental wave length, resonance. The wave pattern is regular with nodes at both boundaries, and the wave amplitude unlimitedly increases with time, as always the case for resonance in linear undamped systems.

Plot3D[vxt[x,t],{x,0,L},{t,0,tMax},PlotPoints->100,AxesLabel->Automatic]

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,t],{x,0,L},{t,0,tMax},PlotPoints->100,AxesLabel->Automatic]

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 ϕ[L, t] == 0. To see it, use Eq. (),  (∂_t)ϕ∝δP. At the open end, the pressure with a good accuracy (although not absolutely) merges with the constant atmospheric pressure, thus δP[L,t]==0. Thus also  (∂_t)ϕ[L,t] == 0. Integrating over time, one obtains ϕ[L, t]==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

(Subscript[λ, 1])=4L

or those of the overtones,

(Subscript[λ, n])=((Subscript[λ, 1])/n),     n=1,3,5,⃛    (even overtones absent)

This is a resonance behavior.

L=10;
λ=4L; (* Wave length.  (Subscript[λ, 1])=4L is the fundamental wave length in this case *)
c=1;     (* Wave speed *)
ω=((2π c)/λ);   (* Frequency *)
tMax=30(L/c);   (* Maximal time of the calculation. L/c is the time for the wave to cross the region *)

Equation=(c^2)(∂_(x,x))ϕ[x,t]-(∂_(t,t))ϕ[x,t]==0;    (* 1d wave equation for the velocity potential *)
IConds={ϕ[x,0]==0,((∂_t)ϕ[x,t]/.t->0)==0};    (* All zero at initial time *)
BConds={  ((∂_x)ϕ[x,t]/.x->0)==Tanh[t]Sin[ ω t], (* The excitation begins gradually to avoid inconsistency of initial and boundary conditions *)
                  ϕ[L,t]==0      (* Closed on the right end *)
      };

Timing[Sol=NDSolve[Join[{Equation},BConds,IConds],ϕ,{x,0,L},{t,0,tMax}, Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid","MinPoints"->100}}];]
ϕxt[x_,t_]:=ϕ[x,t]/.Sol[[1]]
vxt[x_,t_]:=-(∂_xx)ϕxt[xx,t]/.xx->x   (* Velocity of the media in the wave is negative gradient of its potential *)

Calculation with λ=4L - 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,t],{x,0,L},{t,0,tMax},PlotPoints->100,AxesLabel->Automatic]

Further one can check that there is no resonance at λ=2L that would be the second overtone etc. Finding the intensity of the sound emitted from the open end is a more complicated numerical problem because the region external to the pipe is not rectangular. However, there is an analytical solution to this problem.

2d: One-slit diffraction in time

One slit diffraction results in the minima of the wave intensity for the directions specified by the angle θ that satisfy

(Subscript[d, 0]) Sin[θ]=m λ,     m=±1,±2,⃛

Here (Subscript[d, 0]) is the width of the slit and λ is the wave length. We solve the wave equation in a rectangular region with the slip in the middle of the bottom line within the time interval between opening the slit and the wave hitting the far (top) wall. (Stationary wave equation cannot be solved by Mathematica at the moment.)

Lx=10;   Ly=10;   
(Subscript[d, 0])=3;   (* Width of the slit *)
λ=1;    (* Wave length *)
c=1;     (* Wave speed *)
ω=((2π c)/λ);   (* Frequency *)

tMax=10;   (* Has to be chosen so that the first wave front is about to hit the far wall, so that there is no reflection yet. Reflections will complicate the wave pattern *)

(* Modeling the slit *)
PowSlit=50;
FSlit[x_]=(1/(1+(((2x/(Subscript[d, 0])))^(2PowSlit))));
Plot[FSlit[x],{x,-Lx,Lx},PlotRange->All]

(* Creating and solving the equations for the velocity potential *)
(* Boundary conditions: perpendicular component of the velocity on boundaries is zero, except for the slit and is a given periodic function on the slit *)
Equation=(∂_(x,x))ϕ[x,y,t]+(∂_(y,y))ϕ[x,y,t]-(c^(-2))(∂_(t,t))ϕ[x,y,t]==0;
IConds={ϕ[x,y,0]==0,((∂_t)ϕ[x,y,t]/.t->0)==0};    (* All zero at initial time *)
BConds={  ((∂_y)ϕ[x,y,t]/.y->0)==FSlit[x]Tanh[t]Sin[ ω t], (* The slit opens gradually to avoid inconsistency of initial and boundary conditions *)
                 ((∂_y)ϕ[x,y,t]/.y->Ly)==0,(* This BC can be removed *)
                  ((∂_x)ϕ[x,y,t]/.x->Lx)==0,
                   ((∂_x)ϕ[x,y,t]/.x->-Lx)==0      };
Timing[Sol=NDSolve[Join[{Equation},BConds,IConds],ϕ,{x,-Lx,Lx},{y,0,Ly},{t,0,tMax}, Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid","MinPoints"->100}}];]
ϕxyt[x_,y_,t_]:=ϕ[x,y,t]/.Sol[[1]]

tt=tMax;
Plot3D[ϕxyt[x,y,tt],{x,-Lx,Lx},{y,0,Ly},PlotRange->All,AxesLabel->Automatic,PlotPoints->100,AspectRatio->(Ly/(2Lx))]

Plotting the wave pattern together with the theoretical result for the directions of the minima

tt=10;
Show[
DensityPlot[ϕxyt[x,y,tt],{x,-Lx,Lx},{y,0,Ly},PlotRange->All,PlotPoints->100,AspectRatio->(Ly/(2Lx))],
Plot[{Cot[ArcSin[(λ/(Subscript[d, 0]))]]x,-Cot[ArcSin[(λ/(Subscript[d, 0]))]]x,Cot[ArcSin[((2λ)/(Subscript[d, 0]))]]x,-Cot[ArcSin[((2λ)/(Subscript[d, 0]))]]x},{x,-10,10},PlotRange->{{-Ly,Ly},{0,Lx}},AspectRatio->(Lx/(2Ly))]
]

Let us create an animated GIF file to put it onto a web site.

(* Preparing to export into animated GIF *) 
(* Set PlotRange (to avoid jitter), Viewpoint, BoxRatios, and ImageSize *)
tt=tMax;
Plot3D[ϕxyt[x,y,tt],{x,-Lx,Lx},{y,0,Ly},PlotRange->{-0.26,0.26},PlotPoints->100,BoxRatios->{1,(Ly/(2Lx)),0.3},ViewPoint->{0,-100,100},ImageSize->500]

(* Export into animated GIF *)
(* Create a list of frames, Plot3Ds at different times *)

WaveFrame[t_]:=Plot3D[ϕxyt[x,y,t],{x,-Lx,Lx},{y,0,Ly},PlotRange->{-0.26,0.26},PlotPoints->100,BoxRatios->{1,(Ly/(2Lx)),0.3},ViewPoint->{0,-100,100},Axes->None,ImageSize->500]
NFrames=100;
tFrame[k_]:=tMax(k/NFrames);
WaveFrames=Table[WaveFrame[tFrame[k]],{k,1,NFrames}];
Export["D:\\_Lehman\\2010 Fall\\PHY 307\\Notes\Diffraction-one-slit.gif",WaveFrames(*,"TransparentColor"->White*)] ;

2d: Two-slit diffraction in time

Two-slit diffraction results in the maxima of the wave intensity for the directions specified by the angle θ that satisfy

d Sin[θ]=m λ,     m=0,±1,±2,⃛

Here d is the distance between the slits and λ 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.)

Lx=10;   Ly=10;   
d0=1;   (* Width of a slit *)
d=4;     (* Distance between slits *)
λ=1;    (* Wave length *)
c=1;     (* Wave speed *)
ω=((2π c)/λ);   (* Frequency *)

tMax=10;   (* Has to be chosen so that the first wave front is about to hit the far wall, so that there is no reflection yet. Reflections will complicate the wave pattern *)

(* Modeling the slits *)
PowSlit=50;
FSlit[x_]=(1/(1+(((((x+d/2)/(d0/2))))^(2PowSlit))))+(1/(1+(((((x-d/2)/(d0/2))))^(2PowSlit))));
Plot[FSlit[x],{x,-Lx,Lx},PlotRange->All]

(* Creating and solving the equations for the velocity potential *)
(* Boundary conditions: perpendicular component of the velocity on boundaries is zero, except for slits and is a given periodic function on slits *)
Equation=(∂_(x,x))ϕ[x,y,t]+(∂_(y,y))ϕ[x,y,t]-(c^(-2))(∂_(t,t))ϕ[x,y,t]==0;
IConds={ϕ[x,y,0]==0,((∂_t)ϕ[x,y,t]/.t->0)==0};  (* All zero at initial time *)
BConds={  ((∂_y)ϕ[x,y,t]/.y->0)==FSlit[x]Tanh[t]Sin[ ω t], (* Slits open gradually to avoid inconsistency of initial and boundary conditions *)
                  ((∂_y)ϕ[x,y,t]/.y->Ly)==0, (* This BC can be removed *)
                  ((∂_x)ϕ[x,y,t]/.x->Lx)==0,
                   ((∂_x)ϕ[x,y,t]/.x->-Lx)==0      };
Timing[Sol=NDSolve[Join[{Equation},BConds,IConds],ϕ,{x,-Lx,Lx},{y,0,Ly},{t,0,tMax}, Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid","MinPoints"->100}}];]
ϕxyt[x_,y_,t_]:=ϕ[x,y,t]/.Sol[[1]]

tt=10;
Plot3D[ϕxyt[x,y,tt],{x,-Lx,Lx},{y,0,Ly},PlotRange->All,AxesLabel->Automatic,PlotPoints->100,AspectRatio->(Ly/(2Lx))]

Plotting the wave pattern together with the theoretical result for the directions of the maxima

tt=10;
Show[
DensityPlot[ϕxyt[x,y,tt],{x,-Lx,Lx},{y,0,Ly},PlotRange->All,PlotPoints->100,AspectRatio->(Ly/(2Lx))],
Plot[{(10^10)x,Cot[ArcSin[(λ/d)]]x,-Cot[ArcSin[(λ/d)]]x,Cot[ArcSin[((2λ)/d)]]x,-Cot[ArcSin[((2λ)/d)]]x},{x,-10,10},PlotRange->{{-Ly,Ly},{0,Lx}},AspectRatio->(Lx/(2Ly))]
]

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. () falling on a slit from y<0 and solve the time-dependent Schrödinger equation, Eq. () on the other side of the slit in the region 0≤y≤(Subscript[L, y])  and -(Subscript[L, x])≤x≤(Subscript[L, x]). We assume that the region is bounded by rigid walls, so that the boundary condition is Ψ==0 everywhere except the slit of width (Subscript[d, 0]) in the bottom wall, where is Ψ is the same as in the incident plane wave, Ψ~Exp[-iωt]. It is convenient to express the frequency ω through the particle's de Broglie wave length λ as

ω=(ℏ/(2m))(((((2π )/λ)))^2),

see Eqs. () and (). Having λ as a parameter in calculations, we can relate numerical result to the one-slit diffraction formula, Eq. (). We gradually open the slit starting at t=0, to avoid inconsistency of the initial and boundary conditions. As the initial condition we take Ψ[x,y,0] == 0, no particles in the box.

Lx=10;      Ly=10;      (* Sizes of the rectangular spatial region *)
ℏ=1;       (* Planck's constant *)
m=1;      (* Mass of the particle *)
(Subscript[d, 0])=3;    (* Width of the slit *)
λ=1;     (* De Broglie wave length of the particle *)

ω=(ℏ/(2m))(((((2π )/λ)))^2);   (* Frequency corresponding to de Broglie wave length *)

tMax=3;   (* Time of the calculation *)

(* Modeling the slit *)
PowSlit=50;   
FSlit[x_]=(1/(1+(((2x/(Subscript[d, 0])))^(2PowSlit))));
Plot[FSlit[x],{x,-Lx,Lx},PlotRange->All] 
tOpenSlit=0.1;   (* Time to fully open the slit *)

(* Creating and solving the equations for the velocity potential *)
(* Boundary conditions: perpendicular component of the velocity on boundaries is zero, except for the slit and is a given periodic function on the slit *)
Equations={ ℏ(∂_t)ΨRe[x,y,t]+((ℏ^2)/(2m))((∂_(x,x))ΨIm[x,y,t]+(∂_(y,y))ΨIm[x,y,t])==0,
                    -ℏ(∂_t)ΨIm[x,y,t]+((ℏ^2)/(2m))((∂_(x,x))ΨRe[x,y,t]+(∂_(y,y))ΨRe[x,y,t])==0
                    };
IConds={ΨRe[x,y,0]==0,ΨIm[x,y,0]==0};    (* All zero at initial time *)
BConds={  ΨRe[x,0,t]==FSlit[x]Tanh[(t/tOpenSlit)]Cos[ω t],
               ΨIm[x,0,t]==-FSlit[x]Tanh[(t/tOpenSlit)]Sin[ω t],(* The slit opens gradually to avoid inconsistency of initial and boundary conditions *)
                   ΨRe[x,Ly,t]==0,    ΨIm[x,Ly,t]==0,
                  ΨRe[Lx,y,t]==0,     ΨIm[Lx,y,t]==0,
                  ΨRe[-Lx,y,t]==0,    ΨIm[-Lx,y,t]==0  };
Timing[Sol=NDSolve[Join[Equations,BConds,IConds],{ΨRe,ΨIm},{x,-Lx,Lx},{y,0,Ly},{t,0,tMax}, Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid","MinPoints"->80}}];]
ΨRexyt[x_,y_,t_]:=ΨRe[x,y,t]/.Sol[[1]]
ΨImxyt[x_,y_,t_]:=ΨIm[x,y,t]/.Sol[[1]]
Ψ2xyt[x_,y_,t_]:=((ΨRexyt[x,y,t])^2)+((ΨImxyt[x,y,t])^2)

tt=tMax;
Plot3D[Ψ2xyt[x,y,tt],{x,-Lx,Lx},{y,0,Ly},PlotRange->{0,2},PlotPoints->100,BoxRatios->{1,(Ly/(2Lx)),0.3},ViewPoint->{0,-100,100},ImageSize->500]

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 ((|Ψ[r, t]|)^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 ∫((|Ψ[r, t]|)^2)dr 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.

(* Preparing to export into animated GIF *) 
(* Set PlotRange (to avoid jitter), Viewpoint, BoxRatios, and ImageSize *)
WaveFrame[t_]:=Plot3D[Ψ2xyt[x,y,t],{x,-Lx,Lx},{y,0,Ly},PlotRange->{0,2},PlotPoints->100,BoxRatios->{1,(Ly/(2Lx)),0.3},ViewPoint->{0,-100,100},Axes->None,ImageSize->500]
(* Export into animated GIF *)
(* Create a list of frames, Plot3Ds at different times *)
NFrames=100;        tFrame[k_]:=tMax(k/NFrames);
WaveFrames=Table[WaveFrame[tFrame[k]],{k,1,NFrames}];
Export["D:\\_Lehman\\2010 Fall\\PHY 307\\Notes\Diffraction_of_electron_one-slit.gif",WaveFrames(*,"TransparentColor"->White*)] ;