impulse response ode45 help

50 views (last 30 days)
9times6
9times6 on 24 Sep 2022
Edited: David Goodmanson on 30 Sep 2022
How do I implement the dirac delta function while solving an ode involving dirac delta (impulse input)?
Here is my code:
clear all
clc
tspan = 0:0.001:1;
x0=[0 0];
[t1,x1] = ode45(@myfunc,tspan,x0);
plot(t1,x1(:,1))
grid on;
function dxdt=myfunc(t,x)
y = dirac(t)
if y==Inf
y = 1;
end
dx1dt = x(2);
dx2dt = (-500 / 5 )* x(1) + (-10/5)*x(2) +y ;
dxdt = [dx1dt ;dx2dt];
end
It is a simple spring mass damper system. All I get is a flat line in response. I have checked, the Dirac delta function is not working.
  1 Comment
Davide Masiello
Davide Masiello on 24 Sep 2022
Edited: Torsten on 24 Sep 2022
Works for me, I just copied and pasted
clear all
clc
tspan = 0:0.001:1;
x0=[0 0];
[t1,x1] = ode45(@myfunc,tspan,x0);
plot(t1,x1(:,1))
grid on;
function dxdt=myfunc(t,x)
y = dirac(t);
if y==Inf
y = 1;
end
dx1dt = x(2);
dx2dt = (-500 / 5 )* x(1) + (-10/5)*x(2) +y ;
dxdt = [dx1dt ;dx2dt];
end

Sign in to comment.

Answers (1)

Paul
Paul on 24 Sep 2022
Edited: Paul on 24 Sep 2022
Hi 9times6,
Numerical integration schemes like ode45 don't really work for Dirac delta functions.
One option is to compute the response to a step input and then differentiate the result.
tspan = 0:0.001:6;
x0=[0 0];
[t1,x1] = ode45(@myfuncstep,tspan,x0);
x1 = [gradient(x1(:,1),t1) gradient(x1(:,2),t1)];
The more theoretically correct way for this problem is to integrate the differential equation "by hand" from t = 0- to t = 0+. Over this infinitessimal time, the only part of the xdot vector that's integrated is the Dirac delta, and we see that if x(2) is zero at t = 0- then it will integrate to 1 at t = 0+. So we start the integration at t = 0+, but with the non-zero initial condition on x(2), and then we don't have to worry about Dirac input (note that y = 0 in myfunc)
x0=[0 1];
[t2,x2] = ode45(@myfunc,tspan,x0);
Another option is to approximate the Dirac delta with a very narrow rectangular pulse with unit area. In this case we need to force the ode solver to not integrate past the pulse width on the first step.
x0=[0 0];
[t3,x3] = ode45(@myfuncpulse,tspan,x0,odeset('InitialStep',1/1e2));
Finally, for this linear, time invariant system, we can use the Control System Toolbox as a final check.
A = [0 1;-500/5 -10/5];
B = [0;1];
C = [1 0];
[yimp,timp] = impulse(ss(A,B,C,0));
Now compare the results, show they are basically equivalent.
figure
plot(t1,x1(:,1),t2,x2(:,1),t3,x3(:,1),timp,yimp)
grid on;
legend('step+differentiate','initial condition','square pulse','impulse response')
function dxdt=myfuncstep(t,x)
y = 1;
dx1dt = x(2);
dx2dt = (-500 / 5 )* x(1) + (-10/5)*x(2) +y ;
dxdt = [dx1dt ;dx2dt];
end
function dxdt=myfuncpulse(t,x)
y = 1e2*(t <= 1/1e2);
dx1dt = x(2);
dx2dt = (-500 / 5 )* x(1) + (-10/5)*x(2) +y ;
dxdt = [dx1dt ;dx2dt];
end
function dxdt=myfunc(t,x)
y = 0;
dx1dt = x(2);
dx2dt = (-500 / 5 )* x(1) + (-10/5)*x(2) +y ;
dxdt = [dx1dt ;dx2dt];
end
  1 Comment
David Goodmanson
David Goodmanson on 28 Sep 2022
Edited: David Goodmanson on 30 Sep 2022
Hi 9*6
I think the best way to do this is the the second one mentioned by Paul, integrating by hand across the delta function. Then ode45 is provided a nonzero intitial condition for the starting velocity, and the delta function is absent from the myfunc function. However, there is still the scaling of the delta function to consider.
The initial impulse I is the product of a large applied force F0 for a very small duration t0, such that neither of those are known separately, necessarily, but their product I = F0*t0 is known. That's an initial condition you need to specify. Then the initial impulse is I*dirac(t).
The impulse momentum theorem is
F delta_t = I = m delta_v % m = mass
Using F0 and t0 for the left hand side, then (assuming the velocity is zero before the mass gets whacked by the impulse), then
v0 = I/m
is the initial condition v0 for ode45. The use of the dirac function with no multplier and an imposed height of 1 implicitly assumes that I = 1, which is certainly possible, but just one choice among many.

Sign in to comment.

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!