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 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]

share|improve this question
    
You might wish to look into Method -> "EventLocator"... –  J. M. Feb 8 '12 at 18:11
1  
I think you could get more help if you provide relevant parts of your code. –  faleichik Feb 8 '12 at 18:25
    
@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... –  P. Fonseca Feb 8 '12 at 18:33
    
@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... –  P. Fonseca Feb 8 '12 at 18:37
    
@faleichik it correctly detects the first event, but takes no action. The mTotalHead[x] = criticHead[flowrate] + bottomFunction[x] is not being done. Any ideas? –  P. Fonseca Feb 8 '12 at 19:29
show 1 more comment

2 Answers 2

up vote 8 down vote accepted

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

share|improve this answer
add comment

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.

share|improve this answer
    
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. –  P. Fonseca Feb 8 '12 at 19:47
    
Using your example (thank you for it), let's suppose that x[t] cannot be lower than t. –  P. Fonseca Feb 8 '12 at 19:52
    
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. –  David Feb 8 '12 at 19:59
    
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... –  P. Fonseca Feb 8 '12 at 20:06
add comment

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.