I want to model a dynamical system of the form
$\frac{\text{d}x}{\text{d}t} = f(x)+nx\delta(\pi(t-0.2)).
$
The problem is that I have a point source which is reoccurring at fixed time steps (say at 0.2,1.2,2.2…). How can I handle this numerically?
I have tried to figure out solutions and it seems that there are two "easy" approaches:
1) Consider a grid with a fixed step size and define your delta function as 0 on all grid points except the relevant points. However, the influence of the delta function depends strongly on the step size solving this numerically which is why I think this is not correct.
2) Use a Gaussian to represent the Delta-function. This could also work for adaptive step sizes. However, in this case the results depend on the variance of the Gaussian which should be small. Is this a better approach?
If this approach works, I could generate an array consisting of the Gaussians at different time steps which is zero elsewhere and use this with normal ODE-solvers, right?
The third approach is a little bit nasty in the sense that it cannot be used with "normal" ODE solvers. It would be to evaluate f(x) until we approach the point source and take $y(+\epsilon)=e^ny(-\epsilon)$ (Numerical way to deal with Dirac delta.). This does not depend on the variance or step size. Should I go with this one? I could use normal ODE-solvers until the fixed time step and then just take the exponential. And do you have any references on that?
Thank you very much for your help!
Best Answer
If you have times $a_0<a_1<..$ and coefficients $c_k$ and want to solve $$ \dot x(t)=f(x(t))+\sum_{k=0}^N c_kδ(t−a_k) $$ then you can consider the function $u(t)=x(t)-\sum_{k=0}^N c_kH(t−a_k)$ where $H$ is the Heaviside or unit jump/ramp function. Then you get that $$ \dot u(t)=\dot x(t)-\sum_{k=0}^N c_kδ(t−a_k)=f\left(u(t)+\sum_{k=0}^N c_kH(t−a_k)\right) $$ which merely has a discontinuous right side. Numerical solvers of the usual sort with adaptive step size control will not be very happy with that situation, but usually can deal with it in a satisfying manner.
In the modified problem $$ \dot x(t)=f(x(t))+x(t)\sum_{k=0}^N c_kδ(t−a_k) $$ you can again try to combine $\dot x(t)-x(t)\sum_{k=0}^N c_kδ(t−a_k)$ into a derivative of some substitution variable. This is achieved in this case with the help of an integrating factor, resulting in setting $x(t)=\exp\left(\sum_{k=0}^N c_kH(t−a_k)\right)u(t)$ where now $u$ satisfies the ODE $$ \dot u(t)=\exp\left(-\sum_{k=0}^N c_kH(t−a_k)\right)f\left(\exp\left(\sum_{k=0}^N c_kH(t−a_k)\right)u(t)\right). $$ This has again a piecewise continuous right side, so that the solution $u$ is continuous and piecewise continuously differentiable. Again numerical methods can find this solution, in the post-processing the multiplication with the jumpy exponential introduces the jumps for $x$.