Show resonance with wave equation numerically

Question

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.

Answer 1

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}
 ]

Original answer

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.

Answer 2

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}]