Hi there.

I have the following system of ODEs to solve:

$\displaystyle \frac{dx}{dt} = \alpha x^{\frac{3}{2}}$

$\displaystyle \frac{dy}{dt} &=& \sqrt{\frac{\beta}{x^3}}$

Subject to $\displaystyle x(0) = x_0, y(0) = 0 $.

As you can see, the first equation is a first order separable equation which is simple to solve. The second one has $\displaystyle x$ as a variable... so it was my thought to solve the first equation, then subtitute the solution into the 2nd equation and then solve that by separation of variables. So, here's my solution to the first:

$\displaystyle x(t) = \frac{4}{(\alpha t +C)^2}$

Applying the initial condition, I find $\displaystyle C = \pm \frac{2}{\sqrt{x_0}}$.

So I substituted this solution for $\displaystyle x(t)$ into the second differential equation, and set about solving that by separation of variables too. Here's the solution I got:

$\displaystyle y(t) =\frac{\sqrt{\beta}}{32\alpha}\left(\alpha t + C)^4 +D $

Here, $\displaystyle C$ is the integration constant which resulted from the integration of the first differential equation, and $\displaystyle D$ is the integration constant from the integration of the second differential equation. Applying the initial condition for $\displaystyle y$ and substituting in for $\displaystyle C$ as found before, we find:

$\displaystyle D = - \frac{\sqrt{\beta}}{2\alpha x_0^2}$

So, putting that all together, the solved system is:

$\displaystyle x(t) = \frac{4}{(\alpha t \pm \frac{2}{\sqrt{x_0}})^2}$

$\displaystyle y(t) =\frac{\sqrt{\beta}}{32\alpha}\left(\alpha t \pm \frac{2}{\sqrt{x_0}}\right)^4- \frac{\sqrt{\beta}}{2\alpha x_0^2} $

Now, I compared these analytical results against a numerical integration of the same differential equations, and I'm getting confused by the results. Here's the values of the constants and the initial conditions that I used (which are specific to a physical problem I'm trying to solve):

$\displaystyle \alpha = 3.1678\times 10^{-10}$

$\displaystyle \beta = 3.9860 \times 10^5$

$\displaystyle x_0 = 30000 $

$\displaystyle y_0 = 0$

And the integration time interval is [0,1000000].

I integrated with MATLAB's ode45, which is a runge-kutta based integrator, with a maximum step size of 100, and relative and absolute tolerances of $\displaystyle 10^{-6}$.

I've attached the results... 1.png and 2.png show the results of numerical/analytical comparisons for x(t) and y(t) respectively, when I choose $\displaystyle C = + \frac{2}{\sqrt{x_0}}$, and 3.png and 4.png show the same thing, but for $\displaystyle C = -\frac{2}{\sqrt{x_0}}$.

Now, it was my expectation that ONE of the choices of C would produce the correct result, and one would produce the wrong results. However... I have a mix. Because choosing C to be positive gives agreement for y(t), but disagreement for x(t), and visa versa for choosing C to be negative.

How can this be?... surely the constant C must be the same value for both equations, since they describe a single system... how can it be that the first equation requires a negative C, and the second requires a positive C? Have I missed some step in my integration process that should fix this?

Finally, I can't explain why the numerical/analytical match for x(t) is PERFECT (i.e. when correct C is chosen), whereas the match for y(t) diverges with time. Even setting lower integration tolerances and small time steps does not eliminate the error at all, it always diverges at the same rate regardless, leading me to believe that it is not a numerical error, but something inherent in the equations.

Could anybody help shed light on these two issues?