Will Runge-Kutta work for this differential equation?

Question

I am trying to solve this differential equation numerically from $x=0$ to $x=10$.

$$\frac {dy} {dx} = 4 %%%%% y - 5 %%%%%%%% e^{-x }$$

with initial condition $y(0) =1$

I know the exact solution is $y=e^{-x}$.

Can this be solved numerically using the Euler method or the Runge-Kutta method? I have tried it on MATLAB but it doesn't seem to work, but I am wondering if I have made a mistake because apparently RK4 is difficult to break.

Answer 1

The general solution of $y'=4y-5e^{-x}$ is $$ y=Ce^{4x}+e^{-x} $$ which means that once the numerical solution deviates from the exact solution, the error will be magnified with factor $e^{4\Delta x}$.

Why the numerical solution can not stay close to the exact one

There are two sources for the error, the errors from the numerical method and the accumulated floating point noise. Both can be combined in the semi-quantitative error guesstimate $$ (10^{-15}/h+h^4)e^{4 x}. $$ For a fourth order method on a reasonably tame problem one finds thus that $h=\sqrt[5]{\mu}=10^{-3}$ is the optimal step size in double floating point where both error contributions balance. The accumulated error reaches the level $0.01$ of visible deviation in a graph with $y$-range $[-1,1]$ at about $x=\frac14\ln(10^{10})=5.8$. For $h=10^{-4}$ the first term of the floating point accumulation dominates so that the estimated point of visual deviation is even earlier at about $x=\frac14\ln(10^{19})=5.2$.

---

How some accuracy could be achieved over the full interval

To integrate over the full length with at least 3 digits of precision one would need a smaller step size and a number type with a smaller machine constant $\mu$. The error of the optimal step size $h=\sqrt[5]\mu$ is (of the magnitude) $h^4e^{40}=\mu^{4/5}e^{40}$, where $e^{40}=2.35⋅10^{17}$. We want these magnified errors smaller than $10^{−3}$, that is $\mu<(10^{-20})^{5/4}=10^{-25}=2^{-83}$, that is a mantissa length of $25$ decimal or $83$ binary digits.

Thus one would need a quadruple precision or multi-precision data type for a step size $h=10^{−5}$ or smaller for 3-4 valid digits at $x=10$.

---

On the error propagation:

A bound for the local error at some time step $x$ is $$ |f(x)|μ+C(x)h^5 $$ where the first term represents the floating point errors in the evaluation of $f$ and the RK step (assuming no catastrophic cancellation takes place) and the second the error of the method relative to an exact solution close to the current states. $C(x)$ is the coefficient of the first term of the truncation error at $x$ which is an expression in $f$ and its first $4$ derivatives. These errors get summed up to the global error at $x=x_f$ with a magnification factor $e^{L(x_f−x)}$ by the Grönwall lemma, where $L$ is a Lipschitz constant for $f$. This gives approximately (converting one $h$ to $dx$) the integral $$ \int^{x_f}_{x_0}(|f(x)|μ/h+C(x)h^4)e^{L(xf−x)}dx. $$ Then replace the functions in the first term by upper bounds over the region of integration to get an error formula $$ (M_0μ/h+M_4h^4)(e^{L(x_f−x_0)}-1)/L. $$ at $x=x_f$. For tame problems and to get a qualitative picture one can suppose that $M_0/L$ and $M_4/L$ are of magnitude $1$. For more "stiff" problems one has to use more targeted simplifications, considering $Lh$ instead of $h$ is usually a good first step, that is an error guesstimate $(\mu/(Lh)+(Lh)^4)e^{Lx}$ with optimal step size $h=\frac1L\sqrt[5]\mu$.

Answer 2

I have an RK4 code written up, I tested it on there. It seems to work quite well for small values of $x$, but unfortunately explodes after around $6$. See graphs below. Here I used a $h$ value of $1/1000$.

After reducing to $h=1/10000$, I got this, which is slightly better:

Further reductions in $h$ was not possible, and I got a MATLAB error.