3
$\begingroup$

I model the propagation of a pulsating wave. Below is a system of partial differential equations with free boundaries (Neumann conditions), as well as a pulsating component.

Individual components of the field gradient $D[g_x(t, x, y), \{t, 1\}] = D[u(t, x, y), x]$ and $D[g_y(t, x, y), \{t, 1\}] = D[u(t, x, y), y]$ are added for further calculation of the trajectory of a particle moving in this scalar field (there may be an easier way, for now I did it this way).

Below is the code that calculates the solution of the PDE, and then plots the expected trajectory for a single particle that will move in this field.

What is the problem: the field itself is calculated and animated without problems. The problem is in the numerical solution of the ODE system for the generated initial position of the particle in the field and the gradient components. At certain values ​​of time ($t>~7.8$), the solution becomes unstable. A harbinger of this is visible in the plots up to this point - the particle begins to move in the opposite direction. Expected result: the particle will spread along a uniform trajectory further and further from the center of the pulsations.

Below is the animation. Just keep in mind: I chose a lower frame rate, because preparing for the animation takes a significant amount of time, but that's another problem.

PulsatingField dwdt dvdt

(*Define the parameters*) pars = {T = 7}; pulsatingComponent = Sin[3 t] 5 Exp[-5 ((x - 10/2)^2 + (y - 10/2)^2)/5]; (*Define the equations*) eqn = {D[u[t, x, y], {t, 2}] - Laplacian[u[t, x, y], {x, y}] == pulsatingComponent + NeumannValue[-Derivative[1, 0, 0][u][t, x, y], x == 0 || x == 10] + NeumannValue[-Derivative[1, 0, 0][u][t, x, y], y == 0 || y == 10], D[gx[t, x, y], {t, 1}] == D[u[t, x, y], x], D[gy[t, x, y], {t, 1}] == D[u[t, x, y], y]}; (*Initial condition function*) u0[x_, y_] := 5 Exp[-5 ((x - 10/2)^2 + (y - 10/2)^2)/5]; (*Define initial conditions*) ic = {u[0, x, y] == u0[x, y], Derivative[1, 0, 0][u][0, x, y] == 0, gx[0, x, y] == D[u0[x, y], x], gy[0, x, y] == D[u0[x, y], y]}; (*Solve the system of equations*) solutions = NDSolve[{eqn, ic}, {u, gx, gy}, {t, 0, T}, {x, 0, 10}, {y, 0, 10}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}]; (*Calculate the time derivatives for gx and gy*) gxTimeDerivative = D[gx[t, x, y], t] /. solutions; gyTimeDerivative = D[gy[t, x, y], t] /. solutions; (*Generate random initial coordinates for particle*) seed = RandomReal[{1/2, 1/2}, 2]; (*Define the system of ODEs for particle trajectories*) trajectories = NDSolve[{w'[tValue] == First[gxTimeDerivative /. {t -> tValue, x -> w[tValue], y -> v[tValue]}], v'[tValue] == First[gyTimeDerivative /. {t -> tValue, x -> w[tValue], y -> v[tValue]}], w[0] == seed[[1]], v[0] == seed[[2]]}, {w, v}, {tValue, 0, T}]; (*Parametric plot of particle trajectories*) trajectoryPlot = ParametricPlot[ Evaluate[{w[tValue], v[tValue]} /. trajectories], {tValue, 0, T}, PlotRange -> All, AspectRatio -> 1, AxesLabel -> {"x(t)", "y(t)"}, PlotLabel -> "Particle Trajectory Over Time"] (*Plot of w(t) over time*) wPlot = Plot[w[tValue] /. trajectories, {tValue, 0, T}, PlotRange -> All, AxesLabel -> {"t", "w(t)"}, PlotLabel -> "Change of w(t) Over Time", PlotStyle -> Blue] (*Plot of v(t) over time*) vPlot = Plot[v[tValue] /. trajectories, {tValue, 0, T}, PlotRange -> All, AxesLabel -> {"t", "v(t)"}, PlotLabel -> "Change of v(t) Over Time", PlotStyle -> Red]

Code for animation:

(* Compute the values of the derivative w[t,x,y] for different t *) dwdtFrames = Table[gxTimeDerivative /. {t -> tValue}, {tValue, 0, 8, 0.5}]; dvdtFrames = Table[gyTimeDerivative /. {t -> tValue}, {tValue, 0, 8, 0.5}]; (* Create animation for the derivative of w with respect to t *) dwdtAanimatedPlot = ListAnimate[ Table[Plot3D[Evaluate[dwdtFrames[[i]]], {x, 0, 10}, {y, 0, 10}, PlotRange -> All, AxesLabel -> {"x", "y", "\[PartialD]w/\[PartialD]t"}], {i, 1, Length[dwdtFrames]}]]; (* Create animation for the derivative of v with respect to t *) dvdtAanimatedPlot = ListAnimate[ Table[Plot3D[Evaluate[dvdtFrames[[i]]], {x, 0, 10}, {y, 0, 10}, PlotRange -> All, AxesLabel -> {"x", "y", "\[PartialD]v/\[PartialD]t"}], {i, 1, Length[dvdtFrames]}]]; (* Run the animation for the derivative of w *) dwdtAanimatedPlot (* Run the animation for the derivative of v *) dvdtAanimatedPlot
Improve this question
$\endgroup$
2
  • 3
    $\begingroup$ Might be nice to add the animations, to trigger interest. Just a thought. $\endgroup$ Commented Mar 16 at 6:52
  • 1
    $\begingroup$ Solution is good, but "expected result" is wrong. Particles can move in any direction since $\nabla u$ is periodic in time and in space. Please pay attention that your seed = RandomReal[{1/2, 1/2}, 2] is {0.5, 0.5}. $\endgroup$ Commented Mar 17 at 2:22

2 Answers 2

Reset to default
3
$\begingroup$

Well, I copy the solution of the 2d wave equation on a disk with boundary value 0.

Step 1: Implement Bessel function reduction by its 2 step recursion

 BesselReduction = ({Z_[\[Alpha]_?(# >= 2 &),x_] :> (\[Alpha] - 1) 2/x Z[\[Alpha] - 1, x] - Z[\[Alpha] - 2, x] /; MemberQ[ {BesselJ, BesselY, BesselK, BesselI}, Z]});

This is reducing any expression to $J_0,J_1$

 r^5 BesselJ[6, r] -> ( r^5 BesselJ[6, r] //. BesselReduction // FullSimplify) // TraditionalForm

$$r^5 J_6(r)\to 6 \left(3 r^4-128 r^2+640\right) J_1(r)-r \left(r^4-144 r^2+1920\right) J_0(r)$$

Step 2 Define the Helmholtz operator as a function

$$ \text{Helmholtz}\left( \left\{ r \_,\ \phi \_ \right \}, \ k \_, \\" \text{Polar} \\" \right) \ = \ f \longrightarrow -\frac{1}{r} \ \partial_r \left(r\ \partial_r f \right) \ - \ \frac{1}{r^2} \partial_{\phi,\phi} \ f \ - \ k^2\ f $$

Step 3 General particular solution of Helmholtz equation

 DSolve[ Sin[n \[Phi]]^-1 Helmholtz[{r, \[Phi]}, k, "Polar"][ a[r] Sin[n \[Phi]]] == 0, a[r], r][[-1]] /. {C[1] :> 1, C[2] -> 0} {a[r]->BesselJ[n,k r]} 0 == Helmholtz[{r, \[Phi]}, k, "Polar"][ Sin[6 \[Phi]] BesselJ[6, k r]] //. BesselReduction // Simplify True

Step 4 Build the spectral series for boundary condition zero at r=1 for solutions independent of $\phi$

 SpectralSeries = (x/.FindRoot[BesselJ[0,x],{x,#}]&)/@(\[Pi] (Range[120])); SpectralSeries[[1 ;; 5]] {2.40483,5.52008,8.65373,11.7915,14.9309}

Step 5

Define a start function at t=0 and Fourier transform it in the Bessel basis

 Plot[Exp[-(5 r)^2] (1 - r^2)^2, {r, 0, 1}, PlotRange -> All]

start function

 Amp = Quiet[ NIntegrate[r E^-(5 r)^2 (1 - r^2)^2 BesselJ[0, # r], {r, 0, 1}] & /@ SpectralSeries] Amp[[1 ;; 10]] {0.0175069, 0.0139471, 0.00926142, 0.00512205, 0.00235696, 0.000901345,0.000286078, 0.0000752534, 0.0000163832, 2.94784*10^-6}

Step 7

Linear combination of amplitudes, cos of t*Amp^2 and Bessel basis

 f[r_, t_] := Plus @@ ((Amp*(Cos[#^2 t]*BesselJ[0, # r] &) /@ SpectralSeries) RevolutionPlot3D[ f[r, 1.2], {r, 0, 1}, Boxed -> False, Axes -> None ]

Bessel wave angular mode 0

Improve this answer
$\endgroup$
3
$\begingroup$

It's a general problem:

Undamped wave equations propagate even discontiniuites.

By the fundamentals in 1 dimension

 D[u,t,t]-D[u,x,x] == f(t,x)

has the distributional solution

u(r_,x_):= A(t-x) +B(t+x)

as far as A and B have convergent Fourier transforms.

So any numerical error at any time creates an independent new wave.

Summations over space and time diverge after space-time step number * error $\epsilon$ =1 .

This is the reason why only in Fourier space numerical approximations can be sensibly filtered per frequency and wave number range.

So better solve

$$(\omega^2- k^2) u(\omega, k) \ = \ \mathit F(g)(\omega,k))$$ in the space of harmonic waves with apropriate eigenfunctions for the boundary conditions and filter by an appropriate window in $(\omega,k)$ to get a inverse of a finite discrete Fourier sum.

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Can you show some minimal example using Mathematica? $\endgroup$ Commented Mar 16 at 14:26

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.