# Symbolic solution to periodic boundary conditions

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...

• 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$. – Michael Seifert Aug 5 '20 at 21:42
• 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... – Scott G Aug 5 '20 at 22:24

## 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}]] 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]}]; Oy. But remarkably, it does give us a solution that agrees well with the numerical solution: (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[%] 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}}]


However, this solution appears to contain unnatural "shocks" (discontinuities) along the characteristics passing through the "matching point" $$x = 0/1$$. • 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. – Scott G Aug 6 '20 at 16:48
• 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. – Scott G Aug 6 '20 at 16:48
• @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. – Michael Seifert Aug 6 '20 at 17:19
• @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. – Michael Seifert Aug 6 '20 at 17:25