2

I have a second order differential equation $x^2\cdot y''(x)+x\cdot y'(x)+(x^2-10)=0$ with the initial condition $y(100)=1$ and $y'(100)=0$. I want to evaluate $y(x)$ when $x=103$ with 8 digit accuracy.

The first step is rewrite this equation to a system of first order equation. $\begin{cases} y'=w\\ w' = -\frac{w}{x}-\frac{10}{x^2}+1\end{cases}$ with the initial condition $\begin{cases} y(100)=1\\ w(100)=0\end{cases}$. Then I can use fourth-order Runge-Kutta methods to solve the system. In order to achieve 8 digit accuracy, I need to decide the step size $h$. Solving $h^4\le 10^{-8}$, I get $h=0.01$.

Is this seems a reasonable approach? Is there anything special with such a differential equation with the initial condition so that I need use other method to evaluate it?

Qomo
  • 285
  • The error of RK4 is not $h^4$. It is of order $h^4$, meaning that one can bound it by $Ch^4$ where $C$ depends on your equation. I don't have a numerical analysis textbook nearby, so can't give an explicit bound for $C$. –  Mar 03 '13 at 20:31
  • I got hold of a numerical analysis textbook, A first course in numerical analysis by Ralston and Rabinowitz (published by Dover), pages 217-218. The error bounds are a bit too complicated to reproduce here. –  Mar 03 '13 at 22:20
  • @5pm. Thanks for your comment! I looked at the book you mentioned above but still don't get too much insights for estimating C for my example. Do you have any other good references such as examples for solving systems of differential equations without knowing the exact solution? – Qomo Mar 04 '13 at 00:08

2 Answers2

0

It is reasonable to begin with $h=0.01$. Since the theoretical bounds for RK4 are very conservative and also difficult to use, the step halving method may be preferable. It is based on running the computations with two or several step sizes and comparing the results.

For illustration purposes, I used Euler's method (rather than RK4) for this equation. The results were:

  • with $h=0.02$, $y(103)\approx 5.4225$
  • with $h=0.01$, $y(103)\approx 5.4371$
  • with $h=0.005$, $y(103)\approx 5.4443$

Since Euler's method is of 1st order, each result is off by about $Ch$ from the true value. We can estimate $C$ as follows $$ \frac{5.4225-5.4371}{0.02-0.01} \approx -1.46 $$ $$ \frac{5.4371-5.4443}{0.01-0.005} \approx -1.44 $$ The fact that two values of $C$ are close allows us to conclude confidently that the error is bounded by, say, $2h$ when $h$ is small.

For RK4 the error behaves as $Ch^4$, so you would have the differences like $0.02^4-0.01^4$ in the denominator. Try it (and report the results, if you don't mind).

-1

I think that the method of Frobenius is applicable to this problem.

TheWhiteEwok
  • 1