# How do I solve a PDE with a Dirac Delta function?

I have a PDE in the form of $$\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} + u = \delta(x-1),$$ with initial condition $u(x,0)=100$. I'm trying to solve it numerically, but I have no idea on which method should I use. Most of the examples that I refer to has no delta function in the PDE.

Can anyone guide me on what should I do?

Thank you very much.

-
What's the domain? –  Pragabhava Dec 5 '12 at 20:32
Hi there, the domain for x is 0<x<1000. –  Darrel Wee Dec 5 '12 at 20:38
Can we follow the type of the approach in math.stackexchange.com/questions/223552? –  doraemonpaul Dec 6 '12 at 5:01
@doraemonpaul Is exaclty the same thing. See my answer for details. –  Pragabhava Dec 6 '12 at 5:21

You can restate the problem by integrating the equation in $x \in (1-\epsilon,1+\epsilon)$ and taking $\epsilon \to 0$.

$$\int_{1-\epsilon}^{1+\epsilon} \left(\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} + u\right) dx = \int_{1-\epsilon}^{1+\epsilon} \delta(1-x)dx = 1.$$

Now we have to make the assumption that we can exchange the temporal derivative with the integral, and then $$\frac{\partial}{\partial t} \left(\int_{1-\epsilon}^{1+\epsilon} u dx \right) + \int_{1-\epsilon}^{1+\epsilon}\frac{\partial u}{\partial x} dx + \int_{1-\epsilon}^{1+\epsilon} u dx = 1$$ Taking the limit as $\epsilon \to 0$ and assuming it can be exchanged with the time differentiation (and $u$ is integrable I think, I'm not so fresh on my analysis courses) we have that $$\lim_{\epsilon \to 0} \big[u(1+\epsilon, t) - u(1-\epsilon,t)\big] = 1$$ or, in other words, that the jump of the function $u$ at the point $x = 1$ is of magnitude $1$.

Finally, solve for $0 < x < 1$, $1 < x < 1000$ and glue the solutions using the jump condition.

## Analytic solution

For $0 < x < 1$ you have that $$\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} + u = 0$$ with initial conditions $u(x,0) = 100$.

Using the method of characteristics you have the equivalent system \begin{align} \frac{d x}{d \eta} &= 1 &x(\xi)\big|_{\eta = 0} &= \xi\\ \frac{d t}{d \eta} &= 1 &t(\xi)\big|_{\eta = 0} &= 0\\ \frac{d u}{d \eta} &= -u &u(\xi)\big|_{\eta = 0} &= 100 \end{align} wich can easily be solved, leading to \begin{align} x(\xi, \eta) &= \eta + \xi \\ t(\xi,\eta) &= \eta \\ u(\xi,\eta) &= 100e^{-\eta} \end{align} and then we can invert for $(\xi,\eta)$ , leading to the solution $u(x,t) = 100e^{-t}$. On $x=1$, $u(1,t) = 100e^{-t}$, the jump condition states that $$u(1^+,t) = 1 + 100e^{-t}$$ and then,

For $1 < x < 1000$ the problem is $$\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} + u = 0$$ with $u(x,0) = 100$, and $u(1,t) = 1 + 100e^{-t}$. Clearly, there is an inconsistency in $t=0$, $x=1$, derived from $\delta(x-1)$, and the solution is not well defined there. Given that the PDE is a transport equation, the problem will propagate along the characteristic $t = x - 1$, as we will see next.

Case $t < x - 1$.

In this case, the method of characteristics leads to the system \begin{align} \frac{d x}{d \eta} &= 1 &x(\xi)\big|_{\eta = 0} &= \xi\\ \frac{d t}{d \eta} &= 1 &t(\xi)\big|_{\eta = 0} &= \eta\\ \frac{d u}{d \eta} &= -u &u(\xi)\big|_{\eta = 0} &= 100 \end{align} and, as before, the solution is $u(x,t) = 100 e^{-t}$.

Case $x - 1 < t$.

In this case, the method of characteristics leads to the system \begin{align} \frac{d x}{d \xi} &= 1 &x(\eta)\big|_{\xi = 0} &= 1\\ \frac{d t}{d \xi} &= 1 &t(\eta)\big|_{\xi = 0} &= \eta\\ \frac{d u}{d \xi} &= -u &u(\eta)\big|_{\xi = 0} &= 1 + 100e^{-\eta} \end{align} hence \begin{align} x(\xi, \eta) &= \xi + 1 \\ t(\xi,\eta) &= \xi + \eta \\ u(\xi,\eta) &= \left(1 + 100e^{-\eta}\right)e^{-\xi} \end{align} then $\xi = x - 1$, $\eta = t - x +1$ and $$u(x,t) = e^{1-x} + 100e^{-t}$$

It helps to draw a $(x,t)$ diagram to see exactly whats going on. In the whole upper plane, there's the solution $100 e^{-t}$ that comes from the initial condition. Starting from $x = 1$, there is a new component to the solution, namely $e^{1-x}$. This arises from the discontinuity of the function at $x = 1$, and propagates along the characteristics $t = x - c$, where $c \le 1$. Since the solution has to propagate with velocity $1$, the discontinuity wont be seen until time $t = x-1$, where the delta generated wave will reach the observer. The solution can be written as $$u(x,t) = 100e^{-t} + e^{1-x}\big(H(t - x + 1) - H(x - 1)\big)$$

-
Can you elaborate more on how do I solve for 0<x<1 and how to glue the solutions? I have no prior knowledge regarding the jump condition. Thanks. –  Darrel Wee Dec 6 '12 at 8:11
@DarrelWee I've added the complete deduction for the analytic solution. I don't have time to do the numericall algorithm, but the only thing you need to do is, supposing $x_m = 1$, substitute $u_{m+1}$ with $1 - u_m$ in the algorithm. –  Pragabhava Dec 10 '12 at 16:34

Note that this PDE just belongs to the PDE of the form http://eqworld.ipmnet.ru/en/solutions/fpde/fpde1302.pdf so that exact solution can find easily, so you have not necessarily to solve it numerically.

The general solution is $u(x,t)=e^{-x}C(x-t)+e^{-x}\int\delta(x-1)e^x~dx=e^{-x}C(x-t)+e^{1-x}H(x-1)$

$u(x,0)=100$ :

$e^{-x}C(x)+e^{1-x}H(x-1)=100$

$C(x)=100e^x-eH(x-1)$

$\therefore u(x,t)=e^{-x}(100e^{x-t}-eH(x-t-1))+e^{1-x}H(x-1)=100e^{-t}-e^{1-x}H(x-t-1)+e^{1-x}H(x-1)$

-
I wouldn't trust so blindly on such "general" formulas, first and foremost, because contrary to linear ODEs, linear PDEs are tighgtly dependent on initial/boundary conditions (even first order), and second, because that formula most likely has been constructed assuming some regularity on $f$ and $g$. In fact, your solution is wrong, as (using symbolic derivatives) $$u_t + u_x + u = e^{1-x}\delta(x-1).$$ –  Pragabhava Dec 5 '12 at 23:23
This solution looks fine to me. The initial condition is satisfied, and it is straightforward to check that the pde is also satisfied. Application of the product rule of differentiation is valid as the coefficients of the step functions are smooth. Plugging in yields $$u_t+u_x+u=e^{1-x}\delta(x-1)=e^{1-1}\delta(x-1)$$ using the distributional rule $f(x)\delta(x)=f(0)\delta(x)$ for a smooth function $f$. –  user12477 Dec 6 '12 at 0:22
@doraemonpaul's calculation could be dressed up by using a delta net approximation for the delta function (see equations (34)-(40) here), working with $\epsilon>0$ and then showing that the result shown holds when $\epsilon\to 0$, but this doesn't seem necessary to me. –  user12477 Dec 6 '12 at 0:22
@user12477 I didn't knew the distributional rule you mentioned, it might be useful in my reaserch :). Still, I don't believe that just apply this formula is a constructive answer, specially when there are edgy details in the application. The formula relies entirely on the existence and uniqueness of the solution of the equivalent system of ODEs, and it might be clear for you that it can be extended to this case, but it's certainly not trivial. –  Pragabhava Dec 6 '12 at 3:42
Hi there, first and foremost, thanks for all the help. Before this, I have already solved it using the exact solution, which has the same answer as what doraemonpaul did. However, I am trying to validate the answer. Hence, the reason why I tried to use numerical methods to do a comparison between the solutions for the exact and numerical method. Or is there any other better way to validate this solution? Thanks. –  Darrel Wee Dec 6 '12 at 4:13