9
$\begingroup$

I have to integrate an energy along a path. I know the energy at the "beginning" of the path (energy[0]), and I can determine the energy change (gain and loss) along it (energy'[x]). But the values the energy[x] can have are restricted, since, along the path, the energy is always equal or above a certain value (minimumEnergy[x]). This minimum energy limits changes in an abrut manner, and only influence the values downward the path.

If this would be done in a Euler step manner, it would be easy: starting from the beginning of the path, at each new step energy[x]=Max[energy[x-step]+step*energy'[x-step], minimumEnergy[x]].

How can I put this restrictions on NDSolve?

I've tried to artificially change the energy'[x] equation so that it compensates the difference to have the minimumEnergy[x], but this did not went well. Another problem that I have is that energy'[x] is based on the value of energy[x], and the equation has no real value for energy[x] < minimumEnergy[x].

I've tried to figure out how to implement it with StepMonitor, but came to no conclusion.

The system is slightly complex to be placed here, and so I hope the above explanations are good enough for the understanding of my problem.

The following code will not run. I just placed it here for a better understanding. I already implemented J.M. suggestion on the EventLocator.

ans = energy/. Flatten[NDSolve[
    {
        energy[outletChannelEnd] == minimumEnergy[flowrate] + bottomFunction[0],
        energy'[x] == headLoss[flowrate, height4Head[flowrate, energy[x] - bottomFunction[x]]]
    },
    energy,
    {x, 0, 241}, 
    Method -> {
        "EventLocator",
        "Event" -> (energy[x] - (minimumEnergy[flowrate] + bottomFunction[x])),
        "EventAction" :> (energy[x] = minimumEnergy[flowrate] + bottomFunction[x])
    }
]]

The result is as follows, where the blue line is the minimum energy, and the purple is the calculated one:

enter image description here

EDIT: It correctly detects there is an event on the x=73 (I have asked for a print on the event actions), but it doesn't correct the energy'[x]

$\endgroup$
6
  • $\begingroup$ You might wish to look into Method -> "EventLocator"... $\endgroup$ Feb 8, 2012 at 18:11
  • 1
    $\begingroup$ I think you could get more help if you provide relevant parts of your code. $\endgroup$
    – faleichik
    Feb 8, 2012 at 18:25
  • $\begingroup$ @J.M. since there's no definition of energy'[x] when energy[x] < minimumEnergy[x], I get lots of messages from my function (most likely because Mathematica passes this limit when trying to locate the event). But at the end, it looks correct. I'm going to do the simple Euler integration to check results. PS - aren't we glad when we have feedback from the, sometimes, very few messages we add to our own functions... $\endgroup$
    – P. Fonseca
    Feb 8, 2012 at 18:33
  • $\begingroup$ @faleichik I generally do, but I've structure this code so well (I'm very proud of it) that what I have is just calls to bigger functions, and I'm having trouble figuring out how to add it in a leasable way. I'm still thinking... $\endgroup$
    – P. Fonseca
    Feb 8, 2012 at 18:37
  • $\begingroup$ @faleichik it correctly detects the first event, but takes no action. The mTotalHead[x] = criticHead[flowrate] + bottomFunction[x] is not being done. Any ideas? $\endgroup$
    – P. Fonseca
    Feb 8, 2012 at 19:29

2 Answers 2

10
$\begingroup$

You could do something like this which is inspired by the bouncing ball example in the tutorial on the "EventLocator" Method:

g[x_] := Floor[x]

h[x0_, x1_] := Function[{t}, Evaluate@Module[{f},
    Reap[NestWhile[
       Module[{sol, xend},
         sol = First@NDSolve[{f'[x] == 1/2 + g[x] - f[x], 
           f[#] == g[# + .0001]}, f, {x, #, x1},
         Method -> {"EventLocator", "Event" -> (f[x] <= g[x])}];
           xend = sol[[1, 2, 1, 1, -1]];
           Sow[{f[t] /. sol, # < t < xend}];
         xend] &, 
       x0, (# < x1) &],
    _, Piecewise[#2] &][[2, 1]]
]]

The basic idea is to run NDSolve until f[x] drops below some critical function g[x], calculate a new initial value for f, continue the calculation starting from the position where NDSolve stopped until f[x] drops below g[x] again, etc. until a final time x1 is reached. The solution for this example looks like

sol = h[0, 10];
Plot[{sol[x], g[x]}, {x, 0, 10}]

Mathematica graphics

$\endgroup$
0
$\begingroup$

What about piecewise definitions? If your system is of first order and you reach the minimum energy (here: $0.2$), the derivative will vanish for all times identically, therefore this will do the trick:

s=First@NDSolve[
    {
        x'[t] == Piecewise[{{-x[t], x[t] > 0.2}}, 0],
        x[0] == 1
    },
    x,
    {t, 0, 3}
]
Plot[Evaluate[x /. s][t], {t, 0, 3}, PlotRange -> {0, 1}]

enter image description here

By replacing $0.2$ by some time-dependent minimum energy function, you can generalize this according to your needs.

I can only hope that this solves your problem; could you post a minimal (working) example for us to look at otherwise? Your paste just results in a lot of error messages when evaluated.

$\endgroup$
4
  • $\begingroup$ It is not the derivative that has to be restrained. It is the x[t] that can never (will never) be below a certain value, independently of what the derivative says (the certain value is t dependent). So x[t] depends on x’[t], as long as x[t] stays above the minimum value. From that we take that even if x’[t]=0, if the minimum value raises, x[t] will have to raise. $\endgroup$
    – P. Fonseca
    Feb 8, 2012 at 19:47
  • $\begingroup$ Using your example (thank you for it), let's suppose that x[t] cannot be lower than t. $\endgroup$
    – P. Fonseca
    Feb 8, 2012 at 19:52
  • $\begingroup$ I see. Typically it's easier to input mathematical statements into NDSolve, but I'm not sure how to formulate this in your example. Energy conservation is more or less an exterior physical constraint here. $\endgroup$
    – David
    Feb 8, 2012 at 19:59
  • $\begingroup$ The event locator seemed a good idea, but I'm unable to make it work. Each time the x is below the minimumX[t], I would like to reset x[t]=minimumX[t], but it doesn't work, and I don't know why... $\endgroup$
    – P. Fonseca
    Feb 8, 2012 at 20:06

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.