# How to numerically solve a differential equation with a dirac delta function ?

53 views (last 30 days)

Show older comments

Mohit Kumar
on 30 Jun 2020

Commented: Mohit Kumar
on 13 Jul 2020

The differential equation that I want to solve is

Upon using ode45 and the dirac function, the dirac function doesn't seem to have any effect (which makes sense because x never reaches 1 in a numerical solution)

Any ideas on how to solve this numerically?

##### 6 Comments

Alan Stevens
on 1 Jul 2020

Edited: Alan Stevens
on 1 Jul 2020

### Accepted Answer

Walter Roberson
on 1 Jul 2020

Use ode45 to solve the equation over the start time to time 1 (the place the dirac delta is located.) Do not use the if statements like Alan and Mohit show: just end the integration at the point they would take effect. Using if presents theory problems that you can avoid by just stopping at the place of the event.

Now take the final output of that ode45 and give the appropriate kick to the boundary conditions.

Then restart the ode45 from time 1 to the final desired time, passing in the kicked boundary conditions. Do not use if -- again you avoid the theory problem by not having ode45 cross the interval of discontinuity.

##### 11 Comments

David Goodmanson
on 12 Jul 2020

Edited: David Goodmanson
on 12 Jul 2020

Hi Mohit,

there is nothing in the derivation that connects t0 with x0, i.e. it's necessary that x(t0) = x0. Also if two definitie integrals are equal to each other that generally says nothing about the relation of the integrands to each other.

One of the easier ways to show this is by considering the delta function as a tall skinny gaussian function:

delta(x) ~~ (C/abs(a)) *exp(-x^2/a^2) (1)

in the 'limit' as a --> 0. The abs(a) is there because 'a' can be either positive or negative in the a^2 factor but the delta spike is supposed to be positive. (The limit can only be taken after multiplying by some function and integrating, but we will eventually be doing that sometime. It's kind of like Huck Finn saying it's not really stealing if you intend to give it back). Here C = 1/sqrt(pi) is a normalization factor so that the integral of this function = 1.

Consider delta(x(t)-x0)

delta(x(t)-x0) ~~ (C/abs(a))*exp(-(x(t)-x0)^2/a^2)

and suppose x = x0 at t = t0, i.e. x(t0) = x0. This is where the relation between x0 and t0 comes in.

If you expand x(t) in a Taylor series about t=t0 and keep only the constant term and the term that is linear in ((t-t0) you should be able to prove the result. Keep in mind that whatever squared constant appears in the denominator of the exponent, its abs value has to appear in the denominator in the factor in front.

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!