Dirac delta function in maple

Question

I am trying to use a Dirac delta function within a system of equations so that h(t) increases at t=1.5, t=3, t=4.5 etc.

Here is my code below:

A:=30: a:=1: dm:=3: c3:=5: d0:=1/a: t0:=1/dm: h0:=B0: y:=A*a:  cc:=t0*c3:  #same as cc=c3/Bm 
N:=8: T:=0.5:

sys_ode:= diff(h(t),t)=y*sum(Dirac(t-dm*n*T),n=0..N) - exp(1-d(t))*h(t), diff(d(t),t)=exp(1-d(t))*h(t) - cc*d(t);

ics:=h(0)=A*a, d(0)=0:
ND:=dsolve([sys_ode,ics],numeric); #numerical solution to the system


ND(1);

ND(2);

ND(3);

ND(4);

Currently when I run this I get:

ND(1);
      [t = 1., d(t) = HFloat(2.6749858704773097), 

        h(t) = HFloat(23.164506116038023)]
ND(2);
      [t = 2., d(t) = HFloat(2.5365091646635465), 

        h(t) = HFloat(18.95651519442652)]
ND(3);
      [t = 3., d(t) = HFloat(2.376810307084265), 

        h(t) = HFloat(15.018803909414379)]
ND(4);
      [t = 4., d(t) = HFloat(2.1927211646807137), 

        h(t) = HFloat(11.391114874494281)]

But h(t) in theory should be increasing in value since there is input into the system at t=1.5 and at t=3 and not have decreased all the way down to h(t)=11.39 at t=4.

Any ideas where I have gone wrong would be appreciated. Thanks.

Answer 1

Support for numerical integration of ODE problems containing 0-order Dirac functions was added for Maple 2015.0. Since this addition the results from the integration you show look like this:

> ND(1);
       [t = 1., d(t) = 3.01973877409584, h(t) = 37.2561191650856]

> ND(2);
       [t = 2., d(t) = 3.38932165514909, h(t) = 61.6360962313253]

> ND(3);
       [t = 3., d(t) = 3.32743891599543, h(t) = 71.0940887774625]

> ND(4);
       [t = 4., d(t) = 3.59829473587253, h(t) = 79.8444924691185]

The function is indeed increasing.

In prior versions of Maple, no special treatment of Dirac was performed for numeric integration, so unless one hits the Dirac-0 points exactly, they would be ignored, and if they were hit, integration would halt with an undefined. The old way would have worked just as though the Dirac functions were not even there, which is consistent with your result:

> sys_ode:= diff(h(t),t)= - exp(1-d(t))*h(t),
>           diff(d(t),t)=exp(1-d(t))*h(t) - cc*d(t):
> ND:=dsolve([sys_ode,ics],numeric):

> ND(1);
       [t = 1., d(t) = 2.67498587047731, h(t) = 23.1645061160380]

> ND(2);
       [t = 2., d(t) = 2.53650916466355, h(t) = 18.9565151944265]

> ND(3);
       [t = 3., d(t) = 2.37681030708426, h(t) = 15.0188039094144]

> ND(4);
       [t = 4., d(t) = 2.19272116468071, h(t) = 11.3911148744943]