# Confusing Result From System of ODEs...

• Aug 14th 2013, 02:51 PM
Drecks
Confusing Result From System of ODEs...
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?
• Aug 14th 2013, 04:33 PM
HallsofIvy
Re: Confusing Result From System of ODEs...
Quote:

Originally Posted by Drecks
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$

HOW do you get this? Obviously, putting those two different solutions for x into the second equation will give you two different equations for y. Why do you not get two different solutions for y?

Quote:

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?
• Aug 15th 2013, 02:09 AM
Drecks
Re: Confusing Result From System of ODEs...
Quote:

Originally Posted by HallsofIvy
HOW do you get this? Obviously, putting those two different solutions for x into the second equation will give you two different equations for y. Why do you not get two different solutions for y?

There are only two solutions for $\displaystyle x(t)$ due to the fact that $\displaystyle C = \pm 2/\sqrt{x_0}$.

Hence, yes, I do end up with two solutions for $\displaystyle y(t)$ as well, but they only differ due to the $\displaystyle \pm$ value of C as well.

I got to my solution for $\displaystyle y(t)$ by substituting my solution for $\displaystyle x(t)$ into the differential equation for $\displaystyle y(t)$ and then integrating it by separation of variables and substitution. I didn't perform this TWICE with both values of $\displaystyle C$... instead I performed it once with $\displaystyle C$ left as an arbitrary constant, and then considered the two values of $\displaystyle C$ in my solution later on, as follows:

$\displaystyle \frac{dy}{dt} = \sqrt{\frac{\beta}{x(t)^3}}$

$\displaystyle \frac{dy}{dt} = \sqrt{\frac{\beta}{\left(\frac{4}{\left(\alpha t +C\right)^2}\right)^3}}$

$\displaystyle \frac{dy}{dt} = \sqrt{\frac{\beta\left(\alpha t + C\right)^6}{4^3}}$

$\displaystyle \frac{dy}{dt} =\frac{\sqrt{\beta}}{8}\left(\alpha t + C\right)^3$

$\displaystyle dy =\frac{\sqrt{\beta}}{8}\left(\alpha t + C\right)^3 dt$

Make the substitution $\displaystyle u = \alpha t + C$, so $\displaystyle \frac{du}{dt} = \alpha$, so $\displaystyle dt = \frac{du}{\alpha}$.

Hence

$\displaystyle dy = \frac{\sqrt{\beta}}{8} u^3 \frac{du}{\alpha}$

$\displaystyle \int dy = \frac{\sqrt{\beta}}{8\alpha} \int u^3 du$

$\displaystyle y = \frac{\sqrt{\beta}}{8\alpha}u^4 + D$

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

As far as I can tell, nothing has happened to $\displaystyle C$ during that integration process that would change it how it appears in the y(t) solution, compared to how it appears in the x(t) solution. I would expect that choosing $\displaystyle C=+$, would produce either BOTH correct results in both equations, or BOTH wrong results in both equations, and that $\displaystyle C=-$ would do the opposite. I did not expect that $\displaystyle C=+$ would produce the correct result in one equation, and $\displaystyle C=-$ would produce the correct result in the other.
• Aug 15th 2013, 04:09 AM
BobP
Re: Confusing Result From System of ODEs...
I think that you have simply got tangled up with the negative signs.

Separating the variables and integrating the first equation gets you

$\displaystyle -2x^{-1/2}=\alpha t + C,$

so $\displaystyle C$ is going to be negative. No plus or minus.

Then later, when you substitute for $\displaystyle x^{-3/2},$ the substitution will be $\displaystyle x^{-3/2}=-\frac{1}{8}(\alpha t + C)^{3}.$