For an optimization routine I needed to compute the derivative of the right-hand side $\: f_u(x_k, u_k)$ of a discrete-time system $x_{k+1} = f(x_k, u_k)$. Since $\: f_u(x_k, u_k)$ includes terms that are divided by $u_k$, it is not possible to evaluate $f(x_k, u_k=0)$. From the physics, it's pretty clear that a limit must exist. The numerics poses a problem here, though. Any idea, how I could circumvent it?

To illustrate the behavior, I added a double logarithmic plot showing $f_u(x_k, u_k)$ as $u \rightarrow 0$. It doesn't make any sense that $f_u$ escapes here. What the plot does not show is that there are multiple sign changes of $f_u$ as $u$ drops below $10^{-4}$, let's say.

loglog plot of $f_u(x_k, u_k)$ as $u_k \rightarrow 0$

The function $f_u(x_k, u_k)$ is defined as follows:

$$ f_u(x, u) = \partial_u \int_0^T \int_0^t \sin(x_5 + x_6 \tau + u \frac{\tau^2}{2}) \: d\tau. $$

I had a CAS do the job of analytically integrating and derivating the sine term. As already pointed out this produces some division by $u$ which causes the problems when it comes to evaluating at $u=0$.

  • 1
    $\begingroup$ Do you have the analytical form of $f$? If so, can you post it here? $\endgroup$
    – Bill Barth
    Jan 24, 2016 at 15:41
  • 2
    $\begingroup$ I would segregate out the part that involves division by $u$ (which might be smaller than the 2800+ characters long). To that part I would try to apply l'Hopital's Rule to get the limit as $u\to 0$. $\endgroup$
    – hardmath
    Jan 24, 2016 at 17:21
  • 2
    $\begingroup$ I don't think I understand -- is $f(x,u=0)$ finite, or does it go to infinity? In the former case, you can apply l'Hopital's rule, as suggested. In the latter case, however, the function is not continuous at $u=0$, and consequently also not differentiable. $\endgroup$ Jan 25, 2016 at 0:47
  • 4
    $\begingroup$ Can you evaluate $f(x,u)$ for small $u$, or is the error just as bad as for $f_u$? Is the function holomorphic? I feel the question would be much clearer if you said more about the function. Even at 3k characters, can you at least put it into something like a gist (gist.github.com)? $\endgroup$
    – Kirill
    Jan 25, 2016 at 18:00
  • 2
    $\begingroup$ @MaxHerrmann I got the full series for the integrand of $f_u$ in $u$ as $$\sum_{n\geq1} \frac{(-1)^{\lfloor n/2\rfloor}}{2^n (n-1)!}\tau^{2n}Q_n u^{n-1}, $$ where $Q_n$ is $\cos(x_5+x_6\tau)$ (odd $n$), or $\sin(x_5+x_6\tau)$ (even $n$). I expect it's more accurate to take more terms than just the first one. Incidentally, if you've answered your own question with this, please submit a full answer — [SE encourages people to answer their own questions.] $\endgroup$
    – Kirill
    Jan 26, 2016 at 14:40

1 Answer 1


Given that $u \frac{\tau^2}{2} \ll 1$, one way of tackling the numerical oscillations, even before the actual term emerges, is a Taylor approximation in $u$ of the sine term (Thanks to Kirill for the full series):

$$ \sin(\underbrace{x_5 + x_6 \tau}_{u_0} + \underbrace{u \frac{\tau^2}{2}}_{\Delta u}) = \sin(x_5 + x_6 \tau) + \sum_{n=1}^N \frac{(-1)^{\lfloor n/2 \rfloor}}{2^n(n-1)!} \tau^{2n}Q_nu^{n} + o(\Delta u^{N}),$$


$$ Q_n = \begin{cases} \cos(x_5 + x_6 \tau) & n \text{ odd} \\ \sin(x_5 + x_6 \tau) & \text{else.} \end{cases} $$

Introducing this approximation into the original equation and exchanging order of integration and differentiation gives

$$ \begin{array}{ll} f_u(x,u) & = \int_0^T \int_0^t \partial_u \left( \sin(x_5 + x_6 \tau) + \sum_{n=1}^N \frac{(-1)^{\lfloor n/2 \rfloor}}{2^n(n-1)!} \tau^{2n}Q_nu^{n} + o(\Delta u^{N}) \right) \: d\tau \\ & = \int_0^T \int_0^t \sum_{n=1}^N \frac{(-1)^{\lfloor n/2 \rfloor}}{2^n(n-1)!} \tau^{2n}Q_nu^{n-1} + o(\Delta u^{N-1}) \: d\tau. \end{array} $$

A first-order approximation ($N=1$) is then given as

$$ \begin{array}{ll} f_u(x,u) & \approx \int_0^T \int_0^t \cos(x_5 + x_6 \tau) \frac{\tau^2}{2} \: d\tau \\ & = \frac{-3 \cos(x_5)}{x_6^4} + \frac{(6-T^2x_6^2)\cos(x_5 + Tx_6) + 2Tx_6(\sin(x_5) + 2 \sin(x_5 + Tx_6))}{2x_6^4}, \end{array} $$

which does not depend on $u$, naturally.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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