I'm trying to integrate an incredibly simple ODE: $$ y'(x) = -f(y),\quad y(0) = y_0 \ , $$ from $x=0$ to $x=1$. This is a decay type of equation, $f$ is the (variable) decay rate and $y$ is the abundance. Both $y$ and $f(y)$ are positive everywhere.

I'm faced with a couple of constraints. Sometimes $f(y) \gg 1$, which makes this a stiff problem on the timescales I'm interested in. I can evaluate $f(y)$ exactly, but it can be costly. So I'd like to minimize the number of times I have to call it. Luckily, the result doesn't have to be super accurate. This is a tiny routine that will be called often in a much larger code. The current method is obviously horrible: we just take a single forward Euler step: $$ y(1) \approx y(0) + f(y_0) \ . $$ So literally anything suggested will be an improvement :)

Obviously, I could use an RK scheme, probably with adaptive step sizing, or better yet an implicit scheme. However, given the simplicity of the equation I thought there may be something more elegant. For instance, we could rewrite the problem as:

$$ \int_{y_0}^{y_1} \frac{dy}{f(y)} = -1 $$

and attempt to solve for $y_1$, but I'm not sure of the best way to approach this numerically.

Any suggestions appreciated!

  • $\begingroup$ I'm guessing that you don't have any simple way of evaluating $\partial{f}/\partial{y}$? You would need that for an implicit method and since calculation of $f(y)$ is expensive, evaluating using finite differences is unattractive. $\endgroup$ – Bill Greene Apr 30 '15 at 19:06
  • $\begingroup$ Yeah, if it all possible I'd like to avoid relying on $\partial f / \partial y$. I have an explicit expression for $f$, so I could compute the derivative exactly if it turns out to be necessary, but that's a last resort. I feel one should be able to do something explicit with the integral form... $\endgroup$ – Geoff Ryan Apr 30 '15 at 19:38
  • 1
    $\begingroup$ For a scalar ODE there isn't really any point in optimizing things. Just apply a standard (higher-order, adaptive) method. If that isn't fast enough, come back here and let us know. $\endgroup$ – David Ketcheson May 1 '15 at 8:13

Your Answer

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

Browse other questions tagged or ask your own question.