1
$\begingroup$

I wrote a code in Fortran to solve this differential equation using RK4 method:

$$ \frac{dy}{dx}=A\sqrt{\frac{B}{y}+\frac{C}{y^2}} $$

$A$, $B$, and $C$ are some known constants. The problem is that my initial values are $x=0$ and $y=0$ and because $y=0$ is the root of denominator I get a run-time error.

Does anybody know how can I fix this problem?

Improve this question
$\endgroup$
0

3 Answers 3

Reset to default
1
$\begingroup$

Your problem is not well posed. Multiply both sides by $|y|$ to obtain the ODE in the form $$ |y(x)| y'(x) = A \sqrt{By(x)+C}. $$ If your initial condition is $y(0)=0$, then at the initial time you have $$ 0 y'(0) = A\sqrt{C}. $$ This means that if $A\sqrt{C}\neq 0$, no value of $y'(0)$ can satisfy the equation.

In other words, your problem is not that you don't know how to convince whatever ODE integrator you are using to accept your initial condition. Your problem is that your mathematical model makes no sense.

Improve this answer
$\endgroup$
14
  • $\begingroup$ Not necessarily, if at the origin $y \sim \sqrt{x}$ then $y y'$ has a finite limit at the origin. In fact, $y=\sqrt{x}$ is a solution for this ODE, for A=1/2, B=0, and C=1. $\endgroup$
    – Maxim Umansky
    May 18, 2020 at 16:16
  • 2
    $\begingroup$ @MaximUmansky -- true, but is the solution unique? The right hand side is not Lipschitz, and so in general it may not be. As a matter of fact, $y(x)=0$ is also a solution for the same set of coefficients. Non-uniqueness still makes the problem ill-posed as well. $\endgroup$
    – Wolfgang Bangerth
    May 18, 2020 at 19:38
  • $\begingroup$ @ Wolfgang Bangerth Why is $y(x)=0$ a solution here? The left-hand side is zero while the right-hand side is infinity, that does not look right. Perhaps what you mean is $x=0$ for $\forall y$. But that cannot be a solution either because it is not a function. $\endgroup$
    – Maxim Umansky
    May 18, 2020 at 22:15
  • $\begingroup$ Ups, sorry. I meant for $C=0$. For $C=0$, another solution is $y(x)\propto x^p$ with $p=3/2$. $\endgroup$
    – Wolfgang Bangerth
    May 18, 2020 at 22:27
  • $\begingroup$ @ Wolfgang Bangerth Yes, for C=0 there is an analytic solution $x^p$ but $p=2/3$, probably that's what you meant. Still, I don't see non-uniqueness of solutions, with y(0)=0 initial conditions, for given A,B,C. $\endgroup$
    – Maxim Umansky
    May 18, 2020 at 23:23
1
$\begingroup$

The integration curves have the overall form as shown in the figure below, enter image description here they have vertical asymptotes at $y=0$, which certainly creates difficulties for a numerical solver integrating from $y=0$ (x does not matter). This does not necessarily mean the problem is ill posed; that would be the case if the integration curves diverge rapidly as the initial condition is varied, i.e., lack of stability in the sense of Lyapunov. But for the ODE in hand, a simple fix is to step away by a small $\epsilon$ above $y=0$ and integrate from there, the integration curve will stay close to the exact one.

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ A more elengant possibility (in the case the ODE is well posed): use an implicit scheme. In this way, you don't evaluate the function in its singularity. $\endgroup$
    – VoB
    May 18, 2020 at 18:50
  • 2
    $\begingroup$ Yes, an implicit solver would also be an option to avoid the problem with initial conditions. Of course an implicit solver would also have some tolerance $\epsilon$ hidden in it, but this approach does look more elegant. $\endgroup$
    – Maxim Umansky
    May 18, 2020 at 18:56
1
$\begingroup$

If the constants are positive then $y$ is positive and setting $u=y^2$ you get $$u'=2A\sqrt{B\sqrt{|u|}+C}$$ which is a little more regular and allows the initial value $u(0)=0$ without singularity.

This equation is still stiff close to the initial point, which should have consequences for the local error with fixed-step RK4 as long as $u$ is close to zero. One could counter-act this by using smaller step sizes initially.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.