Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am checking a complex situation where a function x1[t1] only exit in a finite time span, and after this a new function x2[t2] will be set up according to the final value of x1, and then repeat the same procedure to build x3[t3], ..., xn[tn].

A most simple 3-nest toy code can be rawly formatted like this

d0 = 1;
{c1, c2, c3} = {1, 2, 3};
end=20;

res1 = NDSolve[{
   x1'[t1] == c1, x1[0] == 0,
   WhenEvent[x1[t1] == d0,
    res2 = NDSolve[{
       x2'[t2] == c2, x2[t1] == x1[t1],
       WhenEvent[x2[t2] == 2*d0,
        res3 = NDSolve[{
           x3'[t3] == c3, x3[t2] == x2[t2]
           }, x3, {t3, t2, end}, MaxSteps -> 1000000];
        "StopIntegration"
        ]}, x2, {t2, t1, end}, MaxSteps -> 1000000];
    "StopIntegration"
    ]}, x1, {t1, 0, end}, MaxSteps -> 1000000];

We can now check res1, res2 and res3, and they work. But I hope that this code can be more concise than its current state, like how we discard the previous functions and code the remaining inner-nest part, especially for multiple-nest cases, ie. a 100-nest case. I think the most proper way would be to go meta-programming together with recursion... am I right?


Update

Based on Albert Retey's answer, we can furthe treat x[t] as a vector function, like this

cVals = {{1, 1/2}, {1/3, 1/4}, {1/5, 1/10}};
tvals = {5, 10, 20};
tstart = 0;
xsol = Quiet[
   NDSolveValue[{
     x'[t] == cVals[[n[t]]], x[tstart] == {0, 0}, n[tstart] == 1,
     WhenEvent[t == tvals[[n[t]]],
      If[TrueQ[n[t] < Length[cVals]], n[t] -> n[t] + 1, "RemoveEvent"]]
     }, x, {t, tstart, 20},
    DiscreteVariables -> {Element[n, Integers]}, MaxSteps -> 1000000
    ], Part::pspec];

Plot[xsol[t], {t, ##}] &[Sequence @@ xsol["Domain"][[1]]]

enter image description here

Where we used cVals = {{1, 1/2}, {1/3, 1/4}, {1/5, 1/10}}.

I think the code can be more general if we can use something like cVals = {1, {1/2, 1/3}, {1/4, 1/5, 1/6}} so that the x[t] can have variable dimensions at different time scopes. I wonder if this idea is feasible?


Update 2

Here is a first try, which shows the idea, but the code does not work yet...

cVals = {1, {1/2, 1/3}};
tvals = {5, 10};
tstart = 0;
xsol = Quiet[
   NDSolveValue[{
     x'[t] == cVals[[n[t]]], x[tstart] == 0, n[tstart] == 1,
     WhenEvent[t > tvals[[n[t]]],
      If[TrueQ[n[t] < Length[cVals]],
       xv = x[t];
       n[t] -> n[t] + 1;
       x[t] -> {xv, xv},
       "RemoveEvent"]]
     }, x, {t, tstart, 20},
    DiscreteVariables -> {Element[n, Integers]}, MaxSteps -> 1000000
    ], Part::pspec];

Plot[xsol[t], {t, ##}] &[Sequence @@ xsol["Domain"][[1]]]
share|improve this question
1  
It's not very hard to write code that does what you want, but why do you think you need to nest these NDSolves? As far as I understand you are just stopping at certain conditions and then run another NDSolve with the final point of the previous as initial conditions, stop that at another condition and so on. No need for a recursion at all. In fact as only a parameter changes in your restarted equations you could handle that within the WhenEvent without even restarting... –  Albert Retey Nov 24 '13 at 19:38
    
@Albert Retey, we can consider the case for this question as a single "path", and the real case I am checking is the coupling among many such paths. During the weekends I have written a code that "mimic" the NDSolve results (as in the question, where the ODE is assumed linear for simplicity) happening at time t1, t2, ... for the multiple paths' coupling. It uses only a single While for dynamic looping. The real problem can be solved by developing an XFEM kink-enrichment method, but for comparison it will be better if we can use finite difference method to solve the problem to some extent. –  saturasl Nov 25 '13 at 20:45
    
@Albert Retey, as you said, it is possible to handle the parameter changing within the WhenEvent without restarting. But how can we do this? Could you provide an example? –  saturasl Nov 25 '13 at 20:54

1 Answer 1

up vote 2 down vote accepted

As mentioned in my comment, it is possible for the example you have shown to achieve the same thing without a recursion and without even restarting NDSolve. The trick is to introduce a discrete variable which is changed at each event. Here is something that I think does the same thing as your code:

cVals = {1, 2, 3};
tstart = 0;
xsol = Quiet[NDSolveValue[{
     x'[t] == cVals[[n[t]]], x[tstart] == 0, n[tstart] == 1,
     WhenEvent[x[t] == n[t]*1,
      If[TrueQ[n[t] < Length[cVals]], n[t] -> n[t] + 1, "RemoveEvent"]
      ]
     },
    x, {t, tstart, 20},
    DiscreteVariables -> {Element[n, Integers]},
    MaxSteps -> 1000000
    ], Part::pspec];

Plot[xsol[t], {t, (xsol@"Domain")[[1, 1]], (xsol@"Domain")[[1, 2]]}]

the code should work for lists of arbitrary length for cVals. The Quiet for Part::pspec is necessary because the equations will be evaluated for symbolic values of n[t]. It could be avoided by defining e.g. cc[n_?NumericQ]:=cVals[[n]] and use that in the first equation.

For other cases you might not be able to reformulate the switch in such a way, but then the first thing I'd try would be a loop, here is code which does that:

c = {1, 2, 3}; d0 = 1; tstart = 0; tend = 20;

solution[n_, nmax_] := Module[{t},
   x[n] = NDSolveValue[{
      x[n]'[t] == c[[n]], x[n][tstop[n - 1]] == x[n - 1][tstop[n - 1]],
      If[n == nmax,
       Unevaluated[Sequence[]],
       WhenEvent[x[n][t] == n*d0, "StopIntegration"]
       ]
      },
     x[n], {t, tstop[n - 1], tend}, MaxSteps -> 1000000
     ];
   tstop[n] = (x[n]@"Domain")[[1, 2]];
   ];

tstop[0] = tstart; x[0][tstart] = 0;

Do[solution[n, Length[c]], {n, 1, Length[c]}]

This checks that the solutions are really calculated as desired:

DownValues[x] // TableForm

Compared to writing the same thing by nesting this won't run into any iteration limit problems. Anyway, for completeness here is code which does what I think is closest to what was asked for originially, a recursive function:

c = {1, 2, 3};
d0 = 1;
x[0][0] = 0;
tend = 20;

solution[n_, nmax_, tprev_] := Module[{t},
   res[n] = NDSolve[{
       x[n]'[t] == c[[n]], x[n][tprev] == x[n - 1][tprev],
       If[n == nmax,
        Unevaluated[Sequence[]],
        WhenEvent[x[n][t] == n*d0, solution[n + 1, nmax, t]; 
         "StopIntegration"]
        ]
       },
      x[n], {t, tprev, tend}, MaxSteps -> 1000000
      ];
   ];

solution[1, Length[c], 0];

DownValues[res]
share|improve this answer
    
This is beautiful code, great! However, I think my code still has its usage, because in coupling case, at interaction points, not only the previous equation should be discarded, we also need to build a sequence of new equations, not a single one. Thus, no one knows exactly how many equations would evently be built before getting the real code running. This is why I was thinking in the direction of meta programming. But your code has already enlightened me, thanks! –  saturasl Nov 27 '13 at 3:02
    
@saturasl: while it is in principle possible to switch off equations for certains periods of time by changing their coefficients so that they become trivial (e.g. x'[t]==0) using the WhenEvent mechanism I also think that there are cases where things become easier and probably more efficient if one restarts NDSolve. I would usually use a normal loop for such cases, where the start condition are determined from the last state of the previous run. I'll probably post code for that when I find a little spare time... –  Albert Retey Nov 28 '13 at 7:21
    
I agree. The switch can be controlled either from outside or inside 'NDSolve'. There should be no nesting-limit for controlling from outside; whereas from inside NDSolve may run into some limit. Thank you for updating. –  saturasl Nov 29 '13 at 23:36
    
I tested and the result shows that your using of normal loop Do is reasonable. As in nonlinear PDE case, each calling of solution will take longer time and so the difference between using Do and Nest or Fold like loops are neglectable. +1 for your example code! And I do wish to see the extra code you mentioned, if you have enough time... –  saturasl Nov 30 '13 at 20:53

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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