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

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

Posted 9 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, 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, 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 9 years ago
 Just to follow up: a bug has been filed on this.
Posted 9 years ago
 Thank you.
Posted 9 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, b == 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 9 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}; 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} /. d[t] -> DiracDelta[t - 10]; ss1 = NDSolveValue[eqs1, s, {t, 0, 100}]; p1 = Plot[ss1[t], {t, 0, 100}, PlotRange -> All] 
Posted 9 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 9 years ago
 Note that DSolve does give the correct solution: In:= eqs = {s'[t] == d[t] - s[t]/10, s == 0, d[t] == DiracDelta[t - 10]}; DSolveValue[eqs, s[t], t] Out= E^(1 - t/10) HeavisideTheta[-10 + t] and Plot[E^(1 - t/10) HeavisideTheta[-10 + t], {t, 0, 100}, PlotRange -> All] gives I have reported the bug....
Posted 9 years ago
 Yes... Style[\[FreakedSmiley], 100] Posted 9 years ago
 Yikes! Posted 9 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, 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...