I have a couple of first order differential equations whose solutions I would like to approximate numerically in my python app.

MATLAB ODEs solvers have a built-in detection of events. I do not see a way to make scipy solvers properly handle differential equations with Dirac delta functions (e.g. like the one below) so I would like to write a simple solver myself.

I know the locations of problematic points beforehand yet I am uncertain on how to properly split a differential equation on a delta function. Take the following simplified example: $$ \frac{dv}{dt} = -\frac{v}{\tau} + \delta(t-2), \quad v(0)=1.0. $$ How should I proceed? Can I just run e.g. Euler method from 0 to 2 including, then add 1 to the last step result (i.e. for $t=2$), and then start a new integration for the remaining time using the result for $t=2$ as an initial value?

EDIT: As explained in the comments below I can only approximate the solutions numerically (I might not know some of the DEs in advance to solve them analytically). So I am interested in how to split a DE on a delta function (in problematic points) while using methods like Euler or Runge-Kutta.


$$ \frac{dv}{dt} +\frac{v}{\tau} = \delta(t-2), \quad v(0)=1 $$ The solution of the associated homogeneous ODE $\quad \frac{dv}{dt} +\frac{v}{\tau}=0\quad $ is $\quad v=C\:e^{-\frac{t}{\tau}}$

General solution of the inhomogeneous ODE on the form $\quad v=f(t)\:e^{-\frac{t}{\tau}}$ $$\frac{dv}{dt} +\frac{v}{\tau} =\:e^{-\frac{t}{\tau}}\left(f'(t)-\frac{f(t)}{\tau} \right)+\frac{1}{\tau}f(t)\:e^{-\frac{t}{\tau}} = \delta(t-2)$$

$$f'(t)= e^{\frac{t}{\tau}}\delta(t-2)$$ $$f(t)=\int e^{\frac{t}{\tau}}\delta(t-2)dt= e^{\frac{2}{\tau}}H(t-2)+c$$ $H$ is the symbol for the Heaviside step function.

$v=e^{-\frac{t}{\tau}}f(t)= e^{-\frac{t}{\tau}}\left( e^{\frac{2}{\tau}}H(t-2)+c\right)$

$$v(t)= e^{-\frac{t-2}{\tau}}H(t-2)+c\:e^{-\frac{t}{\tau}}$$

With condition $v(0)=1 \quad\to\quad 1= e^{-\frac{0-2}{\tau}}H(-2)+c\:e^{-\frac{0}{\tau}}\quad\to\quad c=1- e^{\frac{2}{\tau}}(0)=1$

$$v(t)= e^{-\frac{t-2}{\tau}}H(t-2)+e^{-\frac{t}{\tau}}$$

| cite | improve this answer | |
  • $\begingroup$ JJacquelin: thank you for your answer. It does not solve my problem yet. I know how to solve these kind of DEs analytically (I would use the Laplace transform though). The problem is that I have more complicated cases (all of the form $\frac{dv}{dt} = f(t) + \sum_{c\; \mathrm{in}\; S} \delta(t-c)$) and some I might not even know in advance to use an exact solution. So I have to approximate them using e.g. Euler or Runge-Kutta method where I could use the information of problematic points to avoid huge error. Thus the question is how to properly split a DE on a delta function. $\endgroup$ – Szpilona Jul 12 '16 at 12:17

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.