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

43 次查看（过去 30 天）

显示更早的评论

Mohit Kumar 2020-6-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/557104-how-to-numerically-solve-a-differential-equation-with-a-dirac-delta-function

评论： Mohit Kumar 2020-7-13

采纳的回答： Walter Roberson

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 个评论
显示 4更早的评论隐藏 4更早的评论

Mohit Kumar 2020-7-1

在 MATLAB Online 中打开

Alan, this might be something that you missed in your previous comment. The dirac delta function assumes a value of infinity when its argument goes to zero.

So this creates a lot of approximation in the coarse version of the dirac delta function : "small value", "dirac delta's infinity" and the ode's step size.

I implemented the following code for the dirac delta function:

if abs(x-1)<0.01
        x_dd=-(x+epsilon*(x_d-0.5*(x_d-abs(x_d))*1e5));
        disp(['triggered ' num2str(t)]);
    else
        x_dd=-(x+epsilon*x_d);
end

And used a maximum step with the ode45 as 0.001

However, the dirac delta function triggers only at specific intervals when the function crosses 1, not all of them. This felt extremely weird. A plot is attached.

Thanks

Alan Stevens 2020-7-1

编辑：Alan Stevens 2020-7-1

Hmm. I assumed you just wanted the dxdt - |dxdt| to kick in when x = 1 (The area under the delta function being unity). I'm not sure what you are after if you truly want it to go to infnity (what do you expect the ode function to do with that?). Indeed, if infinity is what you want why bother multiplying it by anything else?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Walter Roberson 2020-7-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/557104-how-to-numerically-solve-a-differential-equation-with-a-dirac-delta-function#answer_459451

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 个评论
显示 9更早的评论隐藏 9更早的评论

David Goodmanson 2020-7-9

编辑：David Goodmanson 2020-7-11

在 MATLAB Online 中打开

Hi Mohit, MODIFIED

I assume that x and t are reversed on your plot, otherwise you have several x values for the same t.

Without the delta function, the solution for x is

x = A*exp(-t/2)*cos(w*t+theta), with w = sqrt(3)/2.

Taking the dellta function term to the right hand side,

x'' +x' +x   =  -abs(xdot)*delta(x-x0)    xdot <0                    (1a)
             =  0                         xdot >0                    (1b)

(It's best to use x0 here and set x0=1 later). There is a negative impulse if x crosses x0 with negative velocity.

If A is not large enough the function will never cross x = x0, so assume that A is greater than x0, maybe much greater than x0.

The delta function in x = x(t) needs to be re-expressed as a delta function in t. A standard identity for the delta function is

delta(x(t)-x0) = (1/abs(( dx/dt | t=t0 ))) * delta(t-t0)
where x(t0) = x0.
[ Here ( dx/dt | t=t0 ) is dx/dt evaluated at t0, and is a constant ].

The second line above assures that as t sweeps through t0, x sweeps through x0. That equation is used to determine t0.

Now

delta(x-x0) = 1/abs(xdot)) *delta(t-t0)

Plugging that into (1a) cancels the xdots and gives

x'' +x' +x  = -delta(t-t0)                     xdot < 0
            = 0                                xdot > 0

where t0 is such that x(t0) = x0.

Under the usual assumptions, x is continuous and xdot is discontinuous across the delta function. For a delta function term of the form

B*delta(t-t0)

then

xdot(t0+eps) - xdot(t0-eps) = B

Since B = -1,

xdot(t0+eps) = xdot(t0-eps) -1

when crossing the delta function. So you can use event detection for x = x0 and then let xdot --> xdot-1 for the new initial condition.

David Goodmanson 2020-7-12

编辑：David Goodmanson 2020-7-12

在 MATLAB Online 中打开

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.

Mohit Kumar 2020-7-13

Ah right, my proof is definitely fallacious. I got carried away trying to contrive the derivation!

Your derivation is nice and concise! Thanks for all the help! Glad to have learnt (although superficially) a hitherto alien concept of the dirac delta function in the context of differential equations.

请先登录，再进行评论。

类别

Mathematics and Optimization Optimization Toolbox Optimization Results Solver Outputs and Iterative Display

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

产品

MATLAB

版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

11 个评论
显示 9更早的评论隐藏 9更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

11 个评论 显示 9更早的评论隐藏 9更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

11 个评论
显示 9更早的评论隐藏 9更早的评论