# Derivative of an ODE marching algorithm such as Euler or Runge Kutta

I have an ordinary differential equation:

$\quad$ $\frac{dP}{dx}$ = $f(x,P,Q)$

which I need to integrate from $x$ = $0$ to $X$, where $P$=$P_0$ at $x$=$0$.

For the moment let's assume I am using simple Euler to solve it:

$\quad$ $P_{i+1}$ = $P_i$ + $h*f(x_i,P_i,Q)$

From the solution to this I get the values of $P_i$ at each $x_i$, and in particular $P_n$ at $x$=$X$. All this works fine where $Q$ is a constant.

I now need to deal with the case where $Q$ is variable and in particular I need to find $\frac{\partial P_n}{\partial Q}$. Obviously I can do it numerically by perturbing $Q$ by a small amount $q$, rerunning the Euler solution to get $P_n'$ to give $\frac{\partial P_n}{\partial Q}$ $\approx$ $\frac{P_n'-P_n}{q}$, however this will take two Euler passes.

So I wish to do it analytically using the first Euler solution. $f$ is an analytic function and so I have all the required derivatives $\frac{\partial f}{\partial x}$, $\frac{\partial f}{\partial P}$, $\frac{\partial f}{\partial Q}$

I also need $\frac{\partial P_n}{\partial P_0}$ although I guess the solution methodology is the same as for $\frac{\partial P_n}{\partial Q}$

Any help doing this would be appreciated.

PS. I am actually using Runge Kutta for this, but am choosing Euler in this question for simplification as the solution technique is all that is important here.

• I doubt the result of this will be nice. The problem is that the perturbation in $Q$ causes perturbations in the prior values of $P$ which then have to be fed into $\frac{\partial f}{\partial P}$. So consider for example (with $P_0$ fixed) $P_2'=P_1'+hf(x_1,P_1',Q')=P_0+hf(x_0,P_0,Q')+hf(x_0,P_0+hf(x_0,P_0,Q'),Q')$. That third term is very annoying to handle, and it only gets worse as you go to higher $n$.
– Ian
Jan 29, 2017 at 16:42
• Not sure I agree. Let's say there is one Euler step. $\quad$ $P_1$=$P_0+h*f(x_0,P_0,Q)$ $\quad$ so $\frac{\partial P_1}{\partial Q}=h*\frac{\partial f}{\partial Q}(x_0,P_0,Q)$ if 2 steps: $\quad$ $P_2$=$P_1+h*f(x_1,P_1,Q)$ So $\quad$ $\frac{\partial P_2}{\partial Q}$=$\frac{\partial P_1}{\partial Q}+h*\frac{\partial f}{\partial Q}(x_1,P_1,Q)$ $\quad$ $=h*\left([\frac{\partial f}{\partial Q}(x_0,P_0,Q)+\frac{\partial f}{\partial Q}(x_1,P_1,Q)\right)$ So it seems to me: $\quad$ $\frac{\partial P_n}{\partial Q}=h*\sum_{i=0}^{n-1}\frac{\partial f}{\partial Q}(x_i,P_i,Q)$ Jan 29, 2017 at 17:37
• sorry I formatted this on separate lines as in my question but it chose to concatenate it all onto one line, hope it is clear Jan 29, 2017 at 17:39
• Your $P_1$ in the perturbed case is not the same as in the unperturbed case. So there is a df/dp term for that.
– Ian
Jan 29, 2017 at 17:47
• I am not perturbing anything. I would only do that for numerical derivatives. I need analytical derivatives. Hence I am using regular differentiation rules. No perturbation. Jan 29, 2017 at 17:55

Using $$\frac{d}{dx}\frac{d}{dQ}P(x,Q)=\frac{d}{dQ}f(x,P(x,Q),Q)=\frac{∂f}{∂P}(x,P,Q)·\frac{dP}{dQ}+\frac{∂f}{∂Q}(x,P,Q)$$ you get for $V=\frac{dP}{dQ}(x,Q)$ the augmented ODE system \begin{align} P'&=f(x,P,Q)\\ V'&=f_P·V+f_Q \end{align} which you can solve like any other ODE system. It also works if $P_0$ is a compontent of $Q$.
• Thanks, I will pursue this approach for $\frac{\partial P_n}{\partial Q}$ What would the equivalent equations be for $\frac{\partial P_n}{\partial P_0}$ ? Jan 30, 2017 at 8:32
• As $f$ does in general not depend directly on $P_0$, you get $W'=f_P·W$ with $W_0=I$ for $W(x)=\frac{∂P(x)}{∂P_0}$. Jan 30, 2017 at 8:40
• @dacfer In a sense you're asking the wrong question. This calculation works for any parameter vector $Q$ (replacing $\frac{d}{dQ}$ by $\nabla_Q$, of course). Initial conditions are parameters too. (Also it should be stressed again that this does not compute the derivative of the numerical solution $P_n$, so writing $\frac{\partial P_n}{\partial Q}$ for this is not appropriate.)
• @Ian. Maybe I am, but I have a requirement to calculate $\frac{\partial P_n}{\partial Q}$ I know I can set $Q=Q_1$, run the n-step ODE solver and get $P_n$ Then I can set $Q=Q_1+\delta$, run it again to get $P_n'$ and estimate the derivative. However I don't want to run the ODE solver twice, particularly as I already have the analytic derivatives for the function at each step. So if I'm asking the wrong question, then feel free to rephrase it, but I know what I need. Jan 31, 2017 at 13:47
• @LutzL I have implemented your equations above using V and W, and in 3 test cases the analytic solution for both $\frac{\partial P_n}{\partial Q}$ and $\frac{\partial P_n}{\partial P0}$ give very similar results to the numerical derivatives calculated with multiple passes of the ODE solver. I had to ensure that the $V$ and $W$ ODE's do not contribute to the step length control algorithm which must only use the base ODE ie. for $P$. Thanks for your solution! Feb 1, 2017 at 10:33