# Show resonance with wave equation numerically

Assume I have a one-dimensional medium of length $$L$$ and in a ring shape. The wave propagation in the medium is described by the wave equation

$$\dfrac{d^2 u}{dt^2} - v^2 \dfrac{d^2 u}{dx^2}=0$$,

with $$v$$ the wave velocity in the medium and the periodic boundary condition $$u(x=0,t) = u(x=L,t)$$.

I know that due to the periodic boundary condition only plane waves $$A \sin(\omega t - kx)$$ with wavenumber $$k = 2\pi m/L$$ ($$m$$ is nonzero integer) and thus with frequency $$\omega = v\times 2\pi m/L$$ can propagate through the medium.

My question is how to show this resonance condition numerically in Mathematica.

I'm imagining doing something like this. I use NDSolve in Mathematica with the wave equation and the periodic boundary condition above. Then somehow I need to tell Mathematica that I inject a sine wave $$A \sin(\omega t)$$ to some location in the medium, say at $$x=0$$. I solve the wave equation for a long time enough to reach to the steady state and measure the oscillation at a new location, $$x=L/3$$ for example. I plot the oscillation amplitude versus the frequency $$\omega$$ and I would expect the oscillation amplitude to be largest at $$\omega = v\times 2\pi m/L$$.

But it seems if I modify the boundary condition at $$x=0$$, it would affect the periodic boundary condition.

• The question is unclear in my view, what do you mean by "can propagate through the medium"? Commented Oct 2, 2021 at 2:32
• @xzczd I mean because of the periodic boundary condition the allowed frequency is quantized into some specific values. So the medium only supports (or supports the most strongly) waves oscillating with those quantized frequencies. Commented Oct 3, 2021 at 0:36
• If I understand what you mean, "I solve the wave equation for a long time enough to reach to the steady state and measure the oscillation at a new location" I don't think this is possible for a wave equation without a damping term, because waves never die down in such medium. Commented Oct 3, 2021 at 2:31
• Ah! You're right. Would you suggest me to add a damping term, \alpha u_t(x,t), to the wave equation? Would this help show the resonance effect? Commented Oct 3, 2021 at 7:03

Here is another attempt. May be this is what is meant.

This plots solution for up to some fixed time. It has Manipulate to allow changing some of the parameters to see the resonance.

To play do

I do not see the amplitude being largest at x=L/3 when w=v*2*Pi*m/L. But may be because time was not large enough. I had to make maximum time fixed, otherwise, NDSolve will solve it again and again each time the slider moved and it become slow.

You can modify and play with it as needed. But I see resonance effect there at this frequency.

code

 Manipulate[
Module[{A, L, k, u, x, t, maxTime = 40},
pde = D[u[x, t], {t, 2}] ==
v^2*D[u[x, t], {x, 2}] +
Piecewise[{{A*Sin[forcingFrequency*t], x == 0}, {0, True}}];
L = 2 Pi;
A = 5;
k = 2*Pi*m/L;
(*set up periodic BC*)
bc = {u[0, t] ==
u[L, t], (D[u[x, t], x] /. x -> 0) == (D[u[x, t], x] /. x -> L)};
bc = {u[0, t] == u[L, t]};

(*set up initial conditons BC. let it all be zero*)
ic = {u[x, 0] == 0, (D[u[x, t], t] /. t -> 0) == 0};

(*solve*)
sol = Quiet@NDSolve[{pde, bc, ic}, u, {x, 0, L}, {t, 0, maxTime}];
Grid[{{Plot3D[
Evaluate[u[x, t] /. sol], {x, 0, 2 Pi}, {t, 0, maxTime},
PlotRange -> All, AxesLabel -> {"x", "time", "u"},
ImageSize -> 300],
Plot[Evaluate[u[L/3, t] /. sol], {t, 0, maxTime},
PlotRange -> {Automatic, {-1, 10}}, ImageSize -> 300,
AxesLabel -> {"time", "Amplitude at x=L/3"}]}
}
]
],
{{forcingFrequency, 1, "forcing Frequency w"}, 0, 10, .1,
Appearance -> "Labeled"},
{{m, 2, "m"}, 0, 10, 1, Appearance -> "Labeled"},
{{v, 2.2, "v"}, 0, 10, .1, Appearance -> "Labeled"},
Grid[{{Button["Set at resonance",
forcingFrequency = v*2*m*Pi/(2*Pi)]}}],
TrackedSymbols :> {forcingFrequency, m, v}
]


I do not understand this

I need to tell Mathematica that I inject a sine wave Asin(ωt) to some location in the medium, say at x=0.

Initial conditions should depend on x only. It can't have t in it. In wave PDE, you give BC, and initial conditions. IC's are initial position and initial velocity. Both which can only depend on x, since t=0.

Unless you mean there is an external source present. But your PDE as you wrote it, does not show this.

You can always plot the solution as follows. (No need to use numerical solver for this, as Mathematica can easily solve it analytically)

ClearAll[u, x, t];
pde = D[u[x, t], {t, 2}] == v^2*D[u[x, t], {x, 2}];
L = 2 Pi; A = 1; w=2;v=1;

(*set up periodic BC*)
bc = {u[0, t] == u[L, t], (D[u[x, t], x] /. x -> 0) == (D[u[x, t], x] /. x -> L)};

(*set up initial conditons BC. Give some initial configuration, and zero speed*)
ic = {u[x, 0] == A*Sin[w*x], (D[u[x, t], t] /. t -> 0) == 0};

(*solve*)
sol = u[x, t] /. First@DSolve[{pde, bc, ic}, u[x, t], {x, t}]

Manipulate[
Plot[sol /. t -> time, {x, 0, 2 Pi}, PlotRange -> {Automatic, {-1, 1}}],
{{time, 0, "time"}, 0, 10, .1, Appearance -> "Labeled"},
TrackedSymbols :> {time}
]


You can modify as needed and add more controls for other parameters if you want.

If you meant that $$\sin(\omega t)$$ is a source, then the PDE should include it. You can modify the code to be

ClearAll[u, x, t];
(*wave PDE with source term*)
pde = D[u[x, t], {t, 2}] == v^2*D[u[x, t], {x, 2}] + A*Sin[w*t];
L = 2 Pi; A = 1; k = 2; w = 2; v = 1;

(*set up periodic BC*)
bc = {u[0, t] ==
u[L, t], (D[u[x, t], x] /. x -> 0) == (D[u[x, t], x] /. x -> L)};

(*set up initial conditons BC. Give some initial configuration, and \
zero speed*)
ic = {u[x, 0] == A*Sin[w*x], (D[u[x, t], t] /. t -> 0) == 0};

(*solve*)
sol = NDSolve[{pde, bc, ic}, u, {x, 0, L}, {t, 0, 10}];
Manipulate[
Plot[Evaluate[u[x, t] /. sol] /. t -> time, {x, 0, 2 Pi},
PlotRange -> {Automatic, {-1, 10}}],
{{time, 0, "time"}, 0, 10, .1, Appearance -> "Labeled"},
TrackedSymbols :> {time}
]


Add controls as need to analyze solution. DSovle could not solve the PDE with the source term added. May be in version 13 it can.

• Hi Nasser! Your answer is helpful. But could you comment on my main question as follows? I would like to show numerically the resonant condition: at $\omega = v*2\pi m/L$ the oscillation is strongest (or the wave is allowed to propagate) and when the frequency is away from this value we see less wave propagation. I'm assuming we place a wave source at the position $x=0$. Appreciate your help. Commented Oct 1, 2021 at 10:15

I'd like to extend my comment to an answer.

I'm imagining doing something like this… I solve the wave equation for a long time enough to reach to the steady state and measure the oscillation at a new location…

I don't think resonance can be illustrated in this manner, because waves never die down in a medium described by wave equation without a damping term, thus there won't be a easy-to-observe steady state. Adding a simple damping term $$\alpha \frac{\partial u}{\partial t}$$ to the equation probably won't help, because it dampens wave no matter what its frequency is.

Then how to illustrate resonance with wave equation? To make it clear, let's make a step backward and revisit how resonance is discussed in textbook. Usually, resonance is illustrated with an ODE like:

$$x''(t)+\omega_0^2 x(t)=F \sin (\omega t)$$

When $$\omega\neq \omega_0$$, the solution is steady:

asol = DSolveValue[{x''[t] + ω0^2 x[t] == F Sin[ω t], x[0] == 1,
x'[0] == 0} /. {ω0 -> 1, ω -> 1/2, F -> 1}, x[t], t] // Simplify
(* Cos[t] + 4/3 Sin[t/2] - (2 Sin[t])/3 *)
Plot[asol, {t, 0, 50}, PlotRange -> All]


While when $$\omega= \omega_0$$, the amplitude becomes larger and larger:

asol2 = DSolveValue[{x''[t] + ω0^2 x[t] == F Sin[ω t], x[0] == 1,
x'[0] == 0} /. {ω0 -> 1, ω -> 1, F -> 1}, x[t], t] // Simplify
(* 1/2 (-((-2 + t) Cos[t]) + Sin[t]) *)
Plot[asol2, {t, 0, 50}, PlotRange -> All]


I believe we should illustrate resonance in a similar manner when dealing with wave equation.

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

With[{u = u[t, x]}, eq = D[u, t, t] == D[u, x, x] + Cos[w t] Exp[-2000 (x - 1/2)^2];
ic = {u == 0, D[u, t] == 0} /. t -> 0;
bc = (u /. x -> 0) == (u /. x -> 1)];

tend = 10;
psol = ParametricNDSolveValue[{eq, ic, bc}, u, {t, 0, tend}, {x, 0, 1}, w,
Method -> mol[400, 2]];


Here Cos[w t] Exp[-2000 (x - 1/2)^2] plays the role of driving force near $$x=\frac{1}{2}$$. Let's check the oscillation at $$x=0.8$$. With a frequency far from $$2n\pi$$, the oscillation is steady:

Plot[psol[1][t, 0.8] // Evaluate, {t, 0, tend}]


But for $$\omega=2\pi$$, resonance happens:

Plot[psol[2 Pi][t, 0.8] // Evaluate, {t, 0, tend}]


• That's very nice. But there is still one thing that bugs me, that is, at resonance the wave keeps growing. The realistic scenario relevant to my question is when I measure transmission of wave through a resonant medium: I send an input wave to the medium at a location and measure the output at a different location. I compare the oscillation amplitudes at two locations, which should tell me that the ratio of two oscillation amplitudes is the largest at resonance. If the output wave keeps growing, the my picture breaks down. This is why I mentioned "steady state" in the first place. Commented Oct 4, 2021 at 1:43