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.

For the following initial value problem

BC[b_]:=x[0]==b;
eq[b_]:=NDSolve[{equation==0, BC[b]}, x, {t,0,1000}];
body[t_,b_]:= x[t]/.eq[b];

one could define the derivative of body with regard to t as

Dbody[t_,b_]:=x'[t]/.eq[b];

However I'm trying to find a way to define this for derivatives of arbitrary order with regard to t, and make the definition of Dbody redundant. Simple invocations of ', D and Derivative failed; for example:

Derivative[1][body][1,1]

returns

body'[1,1]
share|improve this question
    
@m_goldberg: The title does not correspond to what I'm asking. I already get derivatives of arbitrary order wrt t if the initial conditions are constant. It is the case of two variables that I haven't been able to find. –  auxsvr Nov 26 '12 at 17:41
    
eq[b_] := NDSolve[{x'[t] + x[t] == 0, x[0] == b}, x, {t, 0, 1}][[1]] and x'[1/2] /. eq[2] works? –  chris Nov 26 '12 at 18:18
    
@auxsvr. Your original title didn't correspond to your real question. I did my best to fix it. If you don't like what I did, re-edit it, but try to make it clearer this time. –  m_goldberg Nov 26 '12 at 18:30
    
@chris, it works, but it is not exactly what I want. Do you have any idea how to wrap this in a function of {t,b}, so that taking derivatives of it wrt t returns the correct result? –  auxsvr Nov 26 '12 at 20:33
    
does this do what you want ? f[t1_, c_] := Module[{eq}, eq[b_] := NDSolve[{x'[t] + x[t] == 0, x[0] == b}, x, {t, 0, 1}][[1]]; x'[t1] /. eq[c]]; f[1/2, 2] –  chris Nov 26 '12 at 20:52

1 Answer 1

One possibility is to use the StateData` framework of NDSolve[] for the purpose. I'll give a sketch on how one might go about it, and leave the encapsulation as a complete routine to you.

Start with an initial state:

$state = First @ NDSolve`ProcessEquations[{\[FormalX]'[t] + \[FormalX][t] == 0,
                                           \[FormalX][0] == 1}, \[FormalX], t];

(I have chosen the function name to be a formal symbol for safety.)

Whenever you need to propagate your solution, you only need to call NDSolve`Iterate[]:

NDSolve`Iterate[$state, 5]

After propagation (one can go either forward or backward, depending on the need), a call to NDSolve`ProcessSolutions[] yields an InterpolatingFunction[] object...

sol = NDSolve`ProcessSolutions[$state]
   {\[FormalX] -> InterpolatingFunction[{{0., 3.}}, <>]}

...which can now be used for evaluation:

{\[FormalX][2], \[FormalX]'''[3]} /. sol
   {0.135335, -0.0526585}

Why go through all that trouble of using state data, you ask? That's because changing initial conditions is a snap with this framework, by virtue of the function NDSolve`Reinitialize[]:

$state = First[NDSolve`Reinitialize[$state, \[FormalX][0] == -3/2]];

from which you can propagate and evaluate once more:

NDSolve`Iterate[$state, 5];
    sol = NDSolve`ProcessSolutions[$state];
{\[FormalX][2/3], \[FormalX]''[5/2]} /. sol
   {-0.770126, -0.122977}

The gist of it then, is to maintain a state of your DE, propagate if needed before evaluating, and use NDSolve`Reinitialize[] whenever you need to change initial conditions.

share|improve this answer

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.