Absorbing Boundary Condition for Complex/Coupled PDE - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-09-23T06:12:08Z https://mathematica.stackexchange.com/feeds/question/191076 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/191076 7 Absorbing Boundary Condition for Complex/Coupled PDE Questino https://mathematica.stackexchange.com/users/58609 2019-02-07T16:47:13Z 2019-02-12T13:44:18Z <p><strong>Context:</strong> This question is relevant to the physical problem of an excited leaky cavity. The boundary condition on the inside of the cavity is Dirichlet. The other end of the cavity is partially blocked by a potential, and beyond that waves are free to escape to infinity. The dynamics are governed by a modified wave equation.</p> <p>First set the <strong>geometry</strong> of the problem, with an example Gaussian perturbation and potential:</p> <pre><code>{xin, xmax, tmax} = {-300, 400, 3000}; {x0, xp, w, wp, v} = {-150, 0, 2, 8, 1}; c = 0.4; g[x_,t_] := Exp[-(((x - x0) - v*t)/w)^2] gp[x_] := Exp[-((x - xp)/wp)^2] </code></pre> <p>And define the <strong>modified wave equation:</strong></p> <pre><code>weqn = D[psi[x, t], {x, 2}] - D[psi[x, t], {t, 2}] - 2 c*I*D[psi[x, t], {t, 1}] + (c^2 - gp[x])* psi[x, t] == -NeumannValue[-Derivative[0, 1][psi][x, t], x == xmax]; </code></pre> <p>Because waves should be purely outgoing outside of the potential, to avoid integrating out to large x we can implement an absorbing boundary condition by putting in the <code>NeumannValue[...]</code> term above for some xmax sufficiently far from the potential.</p> <p><strong>Initial conditions:</strong></p> <pre><code>ic = {psi[x, 0] == g[x, 0], Derivative[0, 1][psi][x, 0] == Derivative[0, 1][g][x, 0]}; </code></pre> <p><strong>Boundary condition for the inside of the cavity:</strong></p> <pre><code>bcd = {psi[xin, t] == 0}; </code></pre> <p><strong>Attempt at a solution:</strong></p> <p>The absorbing boundary condition seems to require a FEM method, so we write</p> <pre><code>sol = psi /. NDSolve[{weqn, ic, bcd}, psi, {x, xin, xmax}, {t, 0, tmax}, Method -&gt; {"MethodOfLines", "SpatialDiscretization" -&gt; {"FiniteElement", "InterpolationOrder" -&gt; {psi -&gt; 2}, "MeshOptions" -&gt; {"MaxCellMeasure" -&gt; 1}}}][[1, 1]]; </code></pre> <p>This works when <code>c = 0</code>, for which weqn is a simple wave equation with a potential. But when <code>c = 0.4</code> the following warnings are generated,</p> <pre><code>NDSolve::femdpop: The NDSolve`FEM`FEMStiffnessElements operator failed. NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended. </code></pre> <p>The problem seems to be due to the third (imaginary) term in the PDE. So a natural fix would be to split the equation into its real and imaginary parts and solve the coupled PDE, so that all terms are real. In other words, we take <code>psi = psiR + I*psiI</code> and thus obtain</p> <pre><code>weqn1 = D[psiR[x, t], {x, 2}] - D[psiR[x, t], {t, 2}] + 2*c*D[psiI[x, t], {t, 1}] + (c^2 - gp[x])* psiR[x, t] == -NeumannValue[-Derivative[0, 1][psiR][x, t], x == xmax]; weqn2 = D[psiI[x, t], {x, 2}] - D[psiI[x, t], {t, 2}] - 2*c*D[psiR[x, t], {t, 1}] + (c^2 - gp[x])* psiI[x, t] == -NeumannValue[-Derivative[0, 1][psiI][x, t], x == xmax]; ic1 = {psiR[x, 0] == g[x, 0], Derivative[0, 1][psiR][x, 0] == Derivative[0, 1][g][x, 0]}; ic2 = {psiI[x, 0] == 0, Derivative[0, 1][psiI][x, 0] == 0}; bcd1 = {psiR[xin, t] == 0}; bcd2 = {psiI[xin, t] == 0}; </code></pre> <p>Now we solve</p> <pre><code>solP = NDSolve[{weqn1, weqn2, ic1, ic2, bcd1, bcd2}, {psiR, psiI}, {x, xin, xmax}, {t, 0, tmax}, Method -&gt; {"MethodOfLines", "SpatialDiscretization" -&gt; {"FiniteElement", "InterpolationOrder" -&gt; {psiR -&gt; 2, psiI -&gt; 2}, "MeshOptions" -&gt; {"MaxCellMeasure" -&gt; 1}}}][]; </code></pre> <p>Which gives no errors. Then we combine to get the full solution</p> <pre><code>solR = psiR /. solP[]; solI = psiI /. solP[]; sol[x_,t_] := solR[x,t] + I*solI[x,t]; </code></pre> <p>The problem is, now the "absorbing" boundary condition at <code>xmax</code> is partially reflecting, as can be seen in the animation</p> <pre><code>Animate[Plot[Abs@sol[x,t], {x, xin, xmax}, PlotRange -&gt; {All, {0, 0.1}}], {t, 0, 800}] </code></pre> <p><strong>The question is</strong>, how should an absorbing boundary condition be properly implemented for this case of coupled PDEs?</p> <hr> <p><strong>Edit:</strong> To clarify, the problem here is for a leaky cavity with the above PDE on the inside of the potential, but with the simple wave equation on the outside.</p> https://mathematica.stackexchange.com/questions/191076/-/191328#191328 1 Answer by Questino for Absorbing Boundary Condition for Complex/Coupled PDE Questino https://mathematica.stackexchange.com/users/58609 2019-02-11T17:51:29Z 2019-02-12T13:44:18Z <p>We can obtain a solution by modifying the above PDE. To do so, we introduce a smoothing function <code>sm[x]</code> which effectively removes the problematic functions outside of the potential. Then, thanks to @xzczd 's suggestions and the links in the comments, we implement the absorbing boundary condition (abc) for the simple wave equation.</p> <p><strong>Geometry</strong>:</p> <pre><code>{xin, xmax, tmax, pts} = {-60, 100, 300, 350}; {x0, xp, w, wp, v} = {50, 0, 2.5, 2, -1}; c = 0.4; g[x_, t_] := Exp[-(((x - x0) - v*t)/w)^2] gp[x_] := Exp[-((x - xp)/wp)^2] sm[x_] := (-Tanh[2 wp x] + 1)/2 </code></pre> <p><strong>Wave equation and boundary conditions:</strong></p> <pre><code>mol[n_Integer] := {"MethodOfLines", "SpatialDiscretization" -&gt; {"TensorProductGrid", "MaxPoints" -&gt; n, "MinPoints" -&gt; n}} weqn = D[psi[x, t], {x, 2}] - D[psi[x, t], {t, 2}] - 2 c*I*D[psi[x, t], {t, 1}] sm[x] + (c^2 sm[x] - gp[x])* psi[x, t] == 0; ic = {psi[x, 0] == g[x, 0], Derivative[0, 1][psi][x, 0] == Derivative[0, 1][g][x, 0]}; bcd = {psi[xin, t] == 0}; abcr = D[psi[x, t], x] + D[psi[x, t], t] == 0 /. x -&gt; xmax; </code></pre> <p><strong>Solution:</strong></p> <pre><code>nsol = NDSolveValue[{weqn, ic, bcd, abcr}, psi, {x, xin, xmax}, {t, 0, tmax}, Method -&gt; mol[pts]]; </code></pre>