impulse response ode45 help

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.

[Davide Masiello]

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

Answer 1

[Paul]

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

[David Goodmanson]

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.