1
|
5671 Views
|
9 Replies
|
5 Total Likes
View groups...
Share
GROUPS:

# NDSolve - Can anyone see what's wrong in my thinking?

Posted 10 years ago
 Here is a simple differential-algebraic equation. I think that since d[t] has a delta at t=10, s[t] should step up at that point. But it steps down. eqs1 = { s'[t] == d[t] - s[t]/10, s[0] == 0, d[t] == DiracDelta[t - 10] }; ss1 = NDSolveValue[eqs1, s, {t, 0, 100}]; p1 = Plot[ss1[t], {t, 0, 100}, PlotRange -> All]  This did not seem right, but to check I made the sign in front of d[t] in the first equation negative: (* note the - sign in front of d[t] *) eqs2 = { s'[t] == -d[t] - s[t]/10, s[0] == 0, d[t] == DiracDelta[t - 10] }; ss2 = NDSolveValue[eqs2, s, {t, 0, 100}]; p2 = Plot[ss2[t], {t, 0, 100}, PlotRange -> All]  But it does not seem to matter!
9 Replies
Sort By:
Posted 10 years ago
 Strangely the following (which crudely mocks up the delta function with UnitStep functions) works as you expect it to (and as you expect your original system to work): eqs1 = {s'[t] == d[t] - s[t]/10, s[0] == 0, d[t] == UnitStep[t - 10] - UnitStep[t - 11]}; ss1 = NDSolveValue[eqs1, s, {t, 0, 100}]; p1 = Plot[ss1[t], {t, 0, 100}, PlotRange -> All] It looks like you found a bug...
Posted 10 years ago
 Yikes!
Posted 10 years ago
 Yes... Style[\[FreakedSmiley], 100] 
Posted 10 years ago
 Note that DSolve does give the correct solution: In[21]:= eqs = {s'[t] == d[t] - s[t]/10, s[0] == 0, d[t] == DiracDelta[t - 10]}; DSolveValue[eqs, s[t], t] Out[22]= E^(1 - t/10) HeavisideTheta[-10 + t] and Plot[E^(1 - t/10) HeavisideTheta[-10 + t], {t, 0, 100}, PlotRange -> All] givesI have reported the bug....
Posted 10 years ago
 Thanks, David. Interesting. And thanks for reporting it. Regarding my Yikes above, I found one of those on my kitchen wall an hour ago. So it may have escaped.
Posted 10 years ago
 I am not familiar with the code internals but I can easily imagine that a Dirac delta will pose difficulties for any numerical software. There are, however, numerous ways to handle discrete events in NDSolve. There is a very thorough discussion of this in the Documentation Center, under tutorial/NDSolveWhenEvents.One thing that, err, bugs me though. Both of the variants below get the desired result. eqs1 = {s'[t] == DiracDelta[t - 10] - s[t]/10, s[0] == 0}; ss1 = NDSolveValue[eqs1, s, {t, 0, 100}]; p1 = Plot[ss1[t], {t, 0, 100}, PlotRange -> All] eqs1 = {s'[t] == d[t] - s[t]/10, s[0] == 0} /. d[t] -> DiracDelta[t - 10]; ss1 = NDSolveValue[eqs1, s, {t, 0, 100}]; p1 = Plot[ss1[t], {t, 0, 100}, PlotRange -> All] 
Posted 10 years ago
 Thanks, Daniel. And that is interesting. Your code produces correct results when the delta function is included in the single equation system. Mine apparently gets the sign wrong when the delta function is split out into a second equation.I do get your point about numerically handling the delta function. The routine must understand when it has integrated over it, and I would think that would be troublesome for convergence. And it would of course be possible to use a delta function where it escaped from under the integral, and became a monster.I also see that WhenEvent provides a perfectly good way of handling this. The only purpose of the delta under the integral is to increment the solution at that point, which can easily be done with WhenEvent. (WhenEvent was a great addition.)For anyone interested, I ran into this just playing around. Suppose one takes a dose of medication every day at 10PM (2200). The dose moves to the bloodstream with a time constant of 1 hour. In the bloodstream, the concentration decays with a time constant of 7 hours. What is the quantity in the stomach and the blood over time? eqs = { (* stomach *) s'[t] == -s[t]/ts, (* bloodstream *) b'[t] == s[t]/ts - b[t]/tb, s[0] == 0, b[0] == 0, (* dose at 10PM *) WhenEvent[Mod[t - 22, 24] == 0, s[t] -> s[t] + d] }; {ss, bb} = NDSolveValue[eqs /. {ts -> 1, tb -> 7, d -> 1}, {s, b}, {t, 0, 120}]; Plot[{ss[t], bb[t]}, {t, 0, 120}, PlotLegends -> {"s[t]", "b[t]"}, FrameLabel -> {"Time in hours", "Quantity or Concentration"}, PlotTheme -> "Scientific"] 
Posted 10 years ago
 Just to follow up: a bug has been filed on this.
Posted 10 years ago
 Thank you.
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.