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 − c^2 ∆φ@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,y f@x, yD 0
is
f@x, yD = F@xD + G@yD,
where F @ x D and G @ y D are arbitrary functions. The solution of the first-order PDE
∂t f@t, xD − v ∂x f@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
∂ρ ∂t
This is the continuity equation in the differential form. If „V is an infinitesimal volume, ∑ t r„V is the rate of change of the amount (corresponding to r) within this volume, whereas div j „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
∇ P
ρ
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 H v “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
∂t^2 δP − c^2 ∆δP. (4)
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 (3) and integrating it by removing the gradient one obtains
∑ f ∑ t
ã
dP r 0
c^2 r 0
dr.
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 + φ 0 D,
where A is the amplitude, w 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 l are related to w and k by
ω =
2 π T
, k =
2 π λ
One can see that w can be different while k adjusts to w as ck = w.
Let us now find the media velocity in a plane wave: v = − ∇ φ = −A k Sin@ωt − k ⋅ r + φ 0 D.
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 selec- tion 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, tD = φ@rD Sin@ωtD,
or similar with Sin@wtD or Exp@ÂwtD, 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
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 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 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 T 0 while the temperature at the right boundary is 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 D and obtain the ordinary differential equation d^2 T dx^2
The solution of this equation is a linear function of x
T = C 1 x + C 2 = T 0 + HT 1 − T 0 L
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
T 1 − 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
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 e 0 is related to units. The solenoidal part of E is created by the magnetic field B changing in time,
∇ E − (6)
∂ 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 m 0 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
B = ∇ A. (7)
Substituting this into the last Maxwell equation, one obtains the equation for the vector potential A
∇ ∇ A = ∇ H ∇ ⋅ AL − ∆ A μ 0 j + μ 0 ε 0 (8)
∂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 − ∇ φ − (^) (9)
∂ 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
c^2
∂^2 A
∂t^2
− ∆ A μ 0 j −
c^2
∂ ∇ φ ∂t
, c =
μ 0 ε 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. (10).
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 (1) with the electric current density j = r v.
ü Schrödinger equation
In quantum mechanics, dynamics of an isolated system is described by the Schrödinger equation for the complex wave function Y@ r , t D that has the form
— (12)
∂t
H
where Ñ is the Plank constant and 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 D the Hamiltonian function (the energy) is given by
H =
p^2 2 m
Replacing p by the quantum-mechanical operator
p^ ˆ^ = −— ∇ , one obtains the Hamilton operator
H^ ˆ^ = −
—^2 ∆
2 m
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 †Y@ r , t D§^2 ª Y@ r , t D Y@ r , t D*^ 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 station- ary states that depend of time monochromatically,
Y@ r , t D = y@ r D ExpB
ÂE
Ñ
t F,
E being the energy of the particle, and satisfy the stationary Schrödinger equation
H (13)
`
y ã Ey.
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 D = U = const Ψ@r, t D = A Exp@− Hωt − k ⋅ rLD. (14)
Plot3D @ Txt @ x, t D , 8 x, 0, L < , 8 t, 0, tMax < , PlotRange → All, AxesLabel → Automatic D
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 D − κ ∂ x,x T @ x, t D 0 BConds = 8 T @ 0, t D Sin @ω t D , H∂ x T @ x, t D ê. x → L L 0 < ICond = T @ x, 0 D 0;
Sol = NDSolve @ Join @ 8 Equation < , BConds, 8 ICond <D , T, 8 x, 0, L < , 8 t, 0, tMax <D Txt @ x_, t_ D : = T @ x, t D ê. Sol @@ 1 DD TH0,1L@x, tD − 0.1 TH2,0L@x, tD 0
9 T@0, tD Sin@tD, TH1,0L@1, tD 0 = 88 T → InterpolatingFunction@ 88 0., 1.<, 8 0., 20.<<, <>D<<
Plot3D @ Txt @ x, t D , 8 x, 0, L < , 8 t, 0, tMax < , PlotRange → All, AxesLabel → Automatic, PlotPoints → 30 D
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, y D = T 0 , T @ Lx , y D = T 1 , and T @ x , 0D = T A x , Ly E = 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 T 0
HLx − xL y ILy − yM x + HLx − xL y ILy − yM
+ T 1
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.
Now plot the final solution
tt = tMax; Plot3D @ Txyt @ x, y, tt D , 8 x, 0, Lx < , 8 y, 0, Ly < , PlotRange → All, AxesLabel → Automatic D
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 condi- tion 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; T 0 = 1; T 1 = 0.2; Equation = ∂ t T @ x, y, t D − κ I∂ x,x T @ x, y, t D + ∂ y,y T @ x, y, t DM 0; BConds = 8 T @ 0, y, t D T 0 Tanh @ 10 t D , T @ Lx, y, t D T 1 Tanh @ 10 t D , T @ x, 0, t D 0, T @ x, Ly, t D 0 < ; ICond = T @ x, y, 0 D 0;
Sol = NDSolve @ Join @8 Equation < , BConds, 8 ICond <D , T, 8 x, 0, Lx < , 8 y, 0, Ly < , 8 t, 0, tMax < , Method → 8 "MethodOfLines", "SpatialDiscretization" → 8 "TensorProductGrid", "MinPoints" → 100 <<D ; Txyt @ x_, y_, t_ D : = T @ x, y, t D ê. Sol @@ 1 DD
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 @ 8 Txyt @ 0.5, 0.5, t D , Txyt @ 0.25, 0.5, t D , Txyt @ 0.5, 0.3, t D< , 8 t, 0, tMax < , PlotRange → All D
0.2 0.4 0.6 0.8 1.
Relaxation is complete.
tt = tMax; Plot3D @ Txyt @ x, y, tt D , 8 x, 0, Lx < , 8 y, 0, Ly < , PlotRange → All, AxesLabel → Automatic D
Let us now plot the heat flux (or the electric field in the case of the electric potential)
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@ T D 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 T 0. Clearly T 0 ` 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 T 0. 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.
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 + H TIgnition − T0 L Tanh B
t tIgnition
F ;
H∗ Equations and solving ∗L EqT = ∂ t T @ x, t D +
Em c
∂ t m @ x, t D κ ∂ x,x T @ x, t D ;
Eqm = ∂ t m @ x, t D − Γ@ T @ x, t DD m @ x, t D ; IniConds = 8 T @ x, 0 D == T0, m @ x, 0 D 1 < ; BConds = 8 T @ 0, t D == TLeft @ t D , T @ L, t D == T0 < ; solution = NDSolve @ Join @8 EqT, Eqm < , IniConds, BConds D , 8 T, m < , 8 x, 0, L < , 8 t, 0, tMax < , Method → 8 "MethodOfLines", "SpatialDiscretization" → 8 "TensorProductGrid", "MinPoints" → 200 <<D mxt @ x_, t_ D : = m @ x, t D ê. First @ solution D Txt @ x_, t_ D : = T @ x, t D ê. First @ solution D
Out[68]= 88 T^ →^ InterpolatingFunction@^88 0., 1000.<,^8 0., 100.<<,^ <>D, m → InterpolatingFunction@ 88 0., 1000.<, 8 0., 100.<<, <>D<<
In[71]:= Plot3D @ mxt @ x, t D ,^^8 t, 0, 0.5 tMax < ,^^8 x, 0, L < , PlotRange → 8 0, 1 < , PlotPoints → 100, AxesLabel → Automatic D
Out[71]=
One can see a deflagration front propagating with a constant speed from the left to right end of the region.
ü Wave equation
Stationary wave equation for monochromatic waves Df + k^2 f ã 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 1 d region between two walls, 0 § x § L. The left wall makes a small-amplitude sinusoidal motion with a frequency w 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 l = 2 pc ê w, 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 D = v wall@ t D. 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 f and them find the media velocity as v @ x , t D = - ∑ x f@ x , t D. 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 λ 1 = 2 L
or those of the overtones,
λn =
λ 1 n
, n = 1, 2, 3, —
This is a resonance behavior.
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 = c^2 ∂ x,x φ@ x, t D − ∂ t,t φ@ x, t D 0; H∗ 1d wave equation for the velocity potential ∗L IConds = 8 φ@ x, 0 D 0, H∂ t φ@ x, t D ê. t → 0 L 0 < ; H∗ All zero at initial time ∗L BConds = 8 H∂ x φ@ x, t D ê. x → 0 L == Tanh @ t D Sin @ ω t D , H∗ The excitation begins gradually to avoid inconsistency of initial and boundary conditions ∗L H∂ x φ@ x, t D ê. x → L L 0 H∗ Closed on the right end ∗L < ;
Timing @ Sol = NDSolve @ Join @8 Equation < , BConds, IConds D , φ , 8 x, 0, L < , 8 t, 0, tMax < , Method → 8 "MethodOfLines", "SpatialDiscretization" → 8 "TensorProductGrid", "MinPoints" → 100 <<D ; D φ xt @ x_, t_ D : = φ@ x, t D ê. Sol @@ 1 DD vxt @ x_, t_ D : = − ∂ xx φ xt @ xx, t D ê. xx → x H∗ Velocity of the media in the wave is negative gradient of its potential ∗L
8 7.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, t D , 8 x, 0, L < , 8 t, 0, tMax < , PlotPoints → 100, AxesLabel → Automatic D
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.
