How to numerically solve a differential equation with a dirac delta function ?
显示 更早的评论
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 个评论
Alan Stevens
2020-6-30
What are your initial conditions and the time over which you want to solve?
Mohit Kumar
2020-6-30
Alan Stevens
2020-6-30
You could use a coarse approximation to the dirac delta function to see it give a kick to dx/dt. Something like:
d = 0;
if abs(x - 1) < small value
d = (v - abs(v))/2;
end
dXdt = [v; -v - x + d];
However, the "small value" probably needs to be quite large (say 10^-1 or 10^-2 ) to see anything!
Any smaller and ode45 is likely to jump across x = 1 without invoking the condition.
Mohit Kumar
2020-6-30
Mohit Kumar
2020-7-1
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 个投票
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 个评论
Mohit Kumar
2020-7-1
Walter Roberson
2020-7-1
Ah, I was thinking 
was related to time, not to position. To do it for position, you would need to set up an ode event; see ballode https://www.mathworks.com/help/matlab/math/ode-event-location.html#bu7wjcg for a good example.
was related to time, not to position. To do it for position, you would need to set up an ode event; see ballode https://www.mathworks.com/help/matlab/math/ode-event-location.html#bu7wjcg for a good example.Because it is an integral equation, 
will be 0 outside of 
to 
and will be 
inside of 
to 
. 
is not a value, it is a distribution. You cannot substitute infinity at the location where x-1=0 (which is to say, at x=1). Dirac delta is like a probability density function taken in the limit -- and knowing the probability of something does not tell you what the value is. The value associated is whatever the dirac delta is being multiplied by.
will be 0 outside of to and will be inside of to . is not a value, it is a distribution. You cannot substitute infinity at the location where x-1=0 (which is to say, at x=1). Dirac delta is like a probability density function taken in the limit -- and knowing the probability of something does not tell you what the value is. The value associated is whatever the dirac delta is being multiplied by.So, when the event tells you have reached x = 1, then you pull out the last t and last vector of boundary conditions. The vector of boundary conditions will probably be in the order [x, 
] so to calculate the kick, you would take
] so to calculate the kick, you would take boundaries = x{end,:};
kick = boundaries(2) - abs(boundaries(2));
Though at the moment I am not sure which component the kick should be applied to.. the second component, I think (not sure at all.)
Note: in theory x could cross 1 more than once, so you should loop the sequence, not just assume that it will only happen once.
Mohit Kumar
2020-7-3
Walter Roberson
2020-7-3
In the context of differential equations, dirac deltas represent impulses at the time(s) that the dirac parameter becomes zeros, with impulse magnitude given by what the diract delta is being multiplied by.
WIth 
being real-valued, then 
is 0 if 
is positive, and is 
when 
is negative. (The value is then divided by 2.) So there is no impulse if x is increasing ( 
positive ) and a kick downwards if x is decreasing.
being real-valued, then is 0 if is positive, and is when is negative. (The value is then divided by 2.) So there is no impulse if x is increasing ( positive ) and a kick downwards if x is decreasing. The trick will be to figure out what it all means with respect to preserving the equation having to be 0 even under the circumstances of the impulse.
Mohit Kumar
2020-7-5
编辑：Mohit Kumar
2020-7-5
Aha! I understand the problem much better now. I think this might be a viable solution. Let me know of your opinion of this and if this is theoretically valid.
At 
if x is decreasing, the differential equation suddenly becomes 
. As , 
.
if x is decreasing, the differential equation suddenly becomes . As , .So I propose that we can start the simulation from this point such that the output of the differential equation is 
, where 
is whatever value it initially was.
, where is whatever value it initially was.I'll try coding this and update here as to what happens.
David Goodmanson
2020-7-9
编辑：David Goodmanson
2020-7-11
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.
Mohit Kumar
2020-7-10
编辑：Mohit Kumar
2020-7-10
David Goodmanson
2020-7-11
Hi Mohit,
This was a good request of yours, because I realized I made a mistake in the derivation. I modified my comment above to reflect that. For the equation in your comment, take
x'' + x' +x = B*delta(t-t0)
(the delta function is in t, not x) and assume the usual solution with
x continuous across the boundary and x' discontinouos across the boundary.
Integrate both sides from t0-eps to t0+eps and evaluate at the limits to obtain
x'(t0+eps) -x'(t0-eps) +x(t0+eps) -x(t0-eps)
+ Integral{t0-eps,t0+eps} x dt = B
The x(t0+eps) -x(t0-eps) terms cancel because x is continous. The integral term is zero because you are integrating across a vanishly small interval. That leaves
x'(t0+eps) -x'(t0-eps) = B
Mohit Kumar
2020-7-11
Thanks a lot for your response!
For the step that you use here,
delta(x(t)-x0) = (1/abs(( dx/dt | t=t0 ))) * delta(t-t0)
I tried to work out a derivation. Is this the correct jutification?
Hence, equating the functions inside the first and last term of the equality,
David Goodmanson
2020-7-12
编辑：David Goodmanson
2020-7-12
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
更多回答（1 个）
Carlos M. Velez S.
2025-7-24
If you want to apply the Dirac delta function in simulation to continuous-time systems, the following code is enough:
function y = delta_dirac(u)
[n,m] = size(u);
if max(n,m) ==1
dt = 1e-6; % Define a small time increment for the delta function
else
dt = u(2) - u(1);
end
y = zeros(n,m);
for i=1:max(m,n)
if u(i) == 0
y(i) = 1/dt;
else
y(i) = 0;
end
end
1 个评论
Walter Roberson
2025-7-24
ode45() is not a continuous time system, so this function is irrelevant to the situation.
类别
在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
