1
$\begingroup$

Is it possible to get an analytical solution to the following problem? If so, how?

$$a w^{(1,0)}(x,t)+b w^{(0,1)}(x,t)=p(x)+g w(x,t)$$ $$I.C.: w(x,0)=0 | B.C.: w(0,t)=w(1,t) | x \in \left[0, 1\right]$$

Clear["Global`*"]

a = v;
b = 1;
g = -Subscript[\[Lambda], i];
p = Subscript[\[Beta], i] Subscript[\[Psi], 0]*Sin[Pi*x];

ic = w[x, 0] == 0;
bc = w[0, t] == w[1, t];

pde = a*D[w[x, t], x] + b*D[w[x, t], t] == g*w[x, t] + p

DSolve[{pde, ic, bc}, w[x, t], {x, t}]

$$\text{DSolve}\left[\left\{\psi _0 \beta _i \sin (\pi x)=\lambda _i w(x,t)+v w^{(1,0)}(x,t)+w^{(0,1)}(x,t),w(x,0)=0,w(0,t)=w(1,t)\right\},w(x,t),\{x,t\}\right]$$


As you can see, in the current formulation, DSolve is unable to solve this problem with the periodic boundary specified.

Perhaps this...? I have tried to provide some guidance to the solution such as shown below and failed. I would expect the solution to have two solutions, 1) $t > x/v$ and 2) $t <= x/v$ representing the propagation of the wave. Hence, I also tried "helping" there as well to no avail.

Element[x, {0, 1}];
Simplify[DSolve[...],{0<=x<=1}];
Assuming[{0<=x<=1}, DSolve[...]];
Simplify[DSolve[...],{0<=x<=1,t>x/v}];
Assuming[{0<=x<=1,t>x/v}, DSolve[...]];

Note 1: I had a previous question related to simple Mathematica things which are now more understood and I have moved beyond that and spent a good deal of time solving some other problems to build experience, but now I am back to this one with the periodic boundary addition.

Note 2: This post has me worried that periodic boundary conditions are not possible to be solved symbolically as they were specifically not called out...

$\endgroup$
  • $\begingroup$ In this case, even if you set $p(x,t) = 0$ and $g = 0$, you still get no solution. You'd think if Mathematica could handle any equations with periodic BCs analytically, it could handle the ones where the solutions is $w(x,t) = 0$. $\endgroup$ – Michael Seifert Aug 5 '20 at 21:42
  • $\begingroup$ Agreed. I was thinking maybe if one was able to tell the solver the constraints of the variables perhaps the solver would find a solution. I'm not sure if anything I tried was actually did that though... Or perhaps there is some other magic method out there... $\endgroup$ – Scott G Aug 5 '20 at 22:24
5
$\begingroup$

Fourier method

We can attempt to find a periodic solution by decomposing the function $p(x)$ according to its Fourier series; solving the PDE for each term in the Fourier series independently; and then summing these solutions. The result isn't pretty, but it does seem to work.

a = v;
b = 1;
g = -Subscript[\[Lambda], i];
c[n_] = Subscript[\[Beta], i] Subscript[\[Psi], 0] FourierCoefficient[Abs[Sin[\[Pi] x]], x, n, FourierParameters -> {1, -2 \[Pi]}]

Note that FourierCoefficients expects the function to be defined over a region symmetric about the origin. This necessitates the use of $|\sin(\pi x)|$ rather than just $\sin(\pi x)$, so that the function is "correct" over the range $[-\frac12,\frac12]$.

ic = w[x, 0] == 0; 

fourierpde[n_] = a*D[w[x, t], x] + b*D[w[x, t], t] == g*w[x, t] + c[n] Exp[2 \[Pi] I n x]

fouriersoln[n_] = FullSimplify[DSolve[{fourierpde[n], ic}, w[x, t], {x, t}]]

enter image description here

Well, that seems nice enough, at least. Now we just need to add these up.

analyticsoln[x_, t_] = Sum[w[x, t] /. fouriersoln[n], {n, -\[Infinity], \[Infinity]}];

enter image description here

Oy. But remarkably, it does give us a solution that agrees well with the numerical solution:

enter image description here

(Left to right: analytic Fourier series solution, numerical solution, difference between solutions. Note the difference in scale for the third graph.)

A few notes on this method:

  • When I plotted the analytic solution, I did have to wrap analyticsolution[x,t] in Re[] to get a better-quality plot. Without this, the plot has some small "gaps" that I believe are due to rounding error leaving analyticsolution[x,t] with a non-negligible imaginary part.

  • These small imaginary parts could presumably be eliminated by taking FourierCosCoefficient and FourierSinCoefficient separately, solving the PDE for both $p(x) = \sin(2 \pi n x)$ and $\cos(2 \pi n x)$, and then summing. However, I would not be surprised if under this method, Mathematica takes much longer to sum the series (or is unable to find a closed form for the solution at all.)

  • I tried applying FullSimplify to analyticsolution[x,t], but it did not return a simplified result in any reasonable amount of time.

  • An approximate solution, which would be easier for Mathematica to plot and manipulate, could be obtained by truncating the final resummation (i.e. drop all terms with $|n|$ above some threshold. In addition, if you change $p(x)$ and Mathematica is unable to resum the resulting infinite series, then an approximate series solution might be the best you can do.

  • If you wish to solve this equation for inhomogeneous initial conditions $w(x,0) = f(x) \neq 0$, I believe that you could to so as follows: Solve the PDE for $p(x) = 0$ and $w(x,0) = f(x)$; and then add the resulting inhomogeneous source-free solution to the solution obtained via the method above. Note that you will need to periodically extend $f(x)$ over the entire real line in the same way to make this work. You could even, in principle, decompose it in a Fourier series and obtain the solution for inhomogeneous ICs in that way.


Original method (flawed)

My original answer follows. However, the resulting solutions were discontinuous across the characteristics passing through the points $x =$ integer.

Mathematica can be coaxed into providing a solution by extending the source function $p(x)$ to a periodic version that covers the entire real line (with $p(x) = p(x-1)$ for all $x$), and then solving the PDE over the entire real line.

Clear["Global`*"]

a = 1;
b = v;
g = -Subscript[\[Lambda], i];
p = Subscript[\[Beta], i] Subscript[\[Psi], 0] Sin[Pi x] ( 2 HeavisideTheta[Sin[Pi x]] - 1)

ic = w[x, 0] == 0;
bc = w[0, t] == w[1, t];

pde = a*D[w[x, t], x] + b*D[w[x, t], t] == g*w[x, t] + p

DSolve[{pde, ic}, w[x, t], {x, t}]
FullSimplify[%]

enter image description here

The function Sin[Pi x] ( 2 HeavisideTheta[Sin[Pi x]] - 1) gives a rectified sine wave, which has the property that $p(x) = p(x-1)$ for all $x$. A similar type of result (not identical in form, but presumably functionally equivalent) can be found by using

p = Subscript[\[Beta], i] Subscript[\[Psi], 0]*Piecewise[{{Sin[Pi x], Sin[Pi x] >=0}, {- Sin[Pi x], Sin[Pi x] < 0}}]

instead.

However, this solution appears to contain unnatural "shocks" (discontinuities) along the characteristics passing through the "matching point" $x = 0/1$.

enter image description here

Looking at the difference between this solution and the numerical solution (calculated above), the difference appears to be constant in the regions between these characteristics. This suggests that there might be some way to salvage this method by more carefully defining the functions (and/or that there might be a subtle error in my method, or in Mathematica).

$\endgroup$
  • $\begingroup$ I have a typo in the original post, a and b should be switched. Just FYI. It doesn't change the general approach you are suggesting which is awesome. Question though, I've massaged it a bit with FullSimplify[soln, {0 <= x <= 1}] and then put in values Simplify[w[x, t] /. soln2 /. x -> 0.5 /. v -> 0.025 /. Subscript[\[Lambda], i] -> 0.0127 /. Subscript[\[Beta], i] -> 0.0006 /. Subscript[\[Psi], 0] -> 1*^10] to plot the result. $\endgroup$ – Scott G Aug 6 '20 at 16:48
  • $\begingroup$ For $w(x=0.5,t)$ the solution yields this [plot](test) which has a discontinuity when the wave reaches the x position (i.e., when the piecewise solution switches). It seems that this solution is not quite right... I switched a and b for this comment but I haven't updated the original post. $\endgroup$ – Scott G Aug 6 '20 at 16:48
  • $\begingroup$ @ScottG: Interesting. It appears that the "analytical" solution I found differs from the numerical solution (obtained using honest-to-goodness periodic BCs) by a piecewise-constant function. The "steps" between the "pieces" occur at along lines of the form $x = n + vt$ for $n \in \mathbb{Z}$ — or in other words the characteristics passing through the "boundary points". I'll have to think more on what's going on here. $\endgroup$ – Michael Seifert Aug 6 '20 at 17:19
  • $\begingroup$ @ScottG: Other notes: (1) Defining $p(x)$ via a Piecewise function yields the same problem. (2) If $p(x) = \sin (2 \pi x)$ (i.e., no piecewise functions required), the analytical solution and the numerical solution agree. $\endgroup$ – Michael Seifert Aug 6 '20 at 17:25
  • 1
    $\begingroup$ @ScottG: I implemented the Fourier series method and it seems to work much better. See my edited answer. $\endgroup$ – Michael Seifert Aug 6 '20 at 18:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.