3
$\begingroup$

I have a numerically solved a differential equation which has an equilibrium value and a starting point, both included as a parameter. I want to plot the difference between the times it takes to come into a certain proximity of the equilibrium value starting from +/- starting parameter. I'll show you the code and hopefully things get more clear:

bsp[c_, d_] := NDSolve[{x'[t] == -x[t] + d, x[0] == c}, x, {t, 0, 30}]
(* defines the differential eqn *)
bspdiff[d_] := Module[{d1 = d, func, func2, diff, par, x, t, t1},
func := x /. First[bsp[par, d1]];
func2 := x /. First[bsp[-par, d1]];
Plot[
    Abs[(t /. FindRoot[func[t] == 1.01*d1, {t, 2}]) - (t1 /. FindRoot[func2[t1] == 0.99*d1,
        {t1, 2}])],
    {par, 1, 5}]]

As you see, I define two functions, one for the positive, one for the negative value of the starting point. Then I use FindRoot to obtain the times at which they are within 1% of the equilibrium value and deduct those two. I want to plot this result as a function of the starting point. What I get upon executing this module is an empty plot and:

FindRoot::nlnum: The function value {-2.02+x\$70015[2.]} is not a list of numbers with dimensions {1} at {t\$70015} = {2.}. >>

ReplaceAll::reps: {FindRoot[func\$70015[t\$70015]==1.01 2,{t\$70015,2}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

FindRoot::nlnum: The function value {-1.98+x\$70015[2.]} is not a list of numbers with dimensions {1} at {t1\$70015} = {2.}. >>

ReplaceAll::reps: {FindRoot[func2\$70015[t1\$70015]==0.99 2,{t1\$70015,2}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

FindRoot::nlnum: The function value {-1.98+x\$70015[2.]} is not a list of numbers with dimensions {1} at {t1\$70015} = {2.}. >>

General::stop: Further output of FindRoot::nlnum will be suppressed during this calculation. >>

ReplaceAll::reps: {FindRoot[func2\$70015[t1$70015]==0.99 2.,{t1\$70015,2.}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

General::stop: Further output of ReplaceAll::reps will be suppressed during this calculation. >>

When I however execute this manually for single values, it works:

func := x /. First[bsp[2, 1]];
func2 := x /. First[bsp[-2, 1]];
Abs[(t /. FindRoot[func[t] == 1.01, {t, 2}]) - (t1 /. 
FindRoot[func2[t1] == 0.99, {t1, 2}])]

yields a numerical output. What I do not understand - where is my mistake/why does it not work? Or is there a better way to do this? I would be very happy about any help and suggestions.

$\endgroup$

1 Answer 1

4
$\begingroup$

Here's how I would do this, using WhenEvent to find the times you're looking for while solving the differential equation and Sow and Reap to extract them.

Clear[bsp]
bsp[c_ /; c > 1, d_] := With[
  {t1 = First @@ Last@Reap@NDSolve[{
      x'[t] == -x[t] + d
      , x[0] == c
      , WhenEvent[x[t] - 1.01 d == 0, {Sow[t], "StopIntegration"}]
      }
    , x
    , {t, 0, 5000}]
  , t2 = First @@ Last@Reap@NDSolve[{
      x'[t] == -x[t] + d
      , x[0] == -c
      , WhenEvent[x[t] - 0.99 d == 0, {Sow[t], "StopIntegration"}]
    }
    , x
    , {t, 0, 5000}]}
, Abs[t1 - t2]
]

Then, just do

Plot[bsp[c, 1], {c, 1.01, 5}]

resulting in

enter image description here

Pre-Version 9

WhenEvent was introduced in MMA Ver. 9. In previous version, we can use the NDSolve Method called "EventLocator". Here is the adapted code:

Needs["DifferentialEquations`NDSolveUtilities`"];

Clear[bsp]
bsp[c_ /; c > 1, d_] := With[{
  t1 = First @@ Last@Reap@NDSolve[{x'[t] == -x[t] + d, x[0] == c}
     , x
     , {t, 0, 10000}
     , Method -> {"EventLocator"
        , "Event" -> (x[t] == 1.01 d)
        , "EventAction" :> (Sow[t]; "StopIntegration")}
   ]
  , t2 = First @@ Last@Reap@NDSolve[{x'[t] == -x[t] + d, x[0] == -c}
     , x
     , {t, 0, 10000}
     , Method -> {"EventLocator"
        , "Event" -> (x[t] == 0.99 d)
        , "EventAction" :> (Sow[t]; "StopIntegration")}
   ]
  , Abs[t1 - t2]
 ]

Analytic version

Of course, in this particular case, there is an analytic solution to the problem.

sols = x[t] /. First@DSolve[{x'[t] == -x[t] + d, x[0] == c}, x[t], t]
times = t /. {Solve[sols == 1.01 d, t], Solve[sols == 0.99 d /. c -> -c, t]} // Flatten // Quiet
Plot[
  Evaluate[{sols, sols /. c -> -c, 1.01 d, 0.99 d} /. {d -> 1, c -> 2}]
  , {t, 1, 10}
  , GridLines -> {times /. {d -> 1, c -> 2}, None}]

results in

enter image description here

Then, we can define

Clear[bsp]
bsp[c_ /; c > 1, d_] = Abs@*Subtract @@ times;

or, pre MMA-10

bsp[c_ /; c > 1, d_] = Abs~Composition~Subtract @@ times;
$\endgroup$
5
  • $\begingroup$ I am very sorry - I did not include that I am working with an old version of Mathematica (8) as it did not seem relevant to me. I am however very thankful for your answer and might just need to obtain a newer version then ... $\endgroup$
    – Inari
    Commented Sep 1, 2015 at 15:11
  • $\begingroup$ @Inari. A similar method exists for older versions! I will write that up, although I have to remember how it works... $\endgroup$
    – march
    Commented Sep 1, 2015 at 15:15
  • $\begingroup$ @Inari. I believe the second answer will work for you. See updated post. You call bsp[c, d] in the same way as before. $\endgroup$
    – march
    Commented Sep 1, 2015 at 15:36
  • $\begingroup$ I cannot thank you enough. It's wonderful and it works perfectly well. $\endgroup$
    – Inari
    Commented Sep 1, 2015 at 15:51
  • $\begingroup$ @Inari. Thank you for the accept! But it is a good idea to wait for at least a day before accepting just in case someone else wants to post an alternative (and possibly better!) answer. When people see a question with an accepted answer, they less likely to post additional answers; questions with unaccepted answers get more answers! That said, I believe that this is the best way to go about solving this problem. $\endgroup$
    – march
    Commented Sep 1, 2015 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.