Runge-Kutta : error while solving a second order differential equation

Question

I am trying to solve a third order non linear differential equation. I have tried to transform it and I've obtained this problem which is a second order problem:

I am trying to implement a fourth order Range-Kutta algorithm in order to solve it by writing it like this :

Here is my code for the Range-Kutta algorithm :

import numpy as np
import matplotlib.pyplot as plt

''''X,Y = integrate(F,x,y,xStop,h).
4th-order Runge-Kutta method for solving the initial value problem {y}' = {F(x,{y})}, where {y} = {y[0],y[1],...,y[n-1]}.
x,y = initial conditions
xStop = terminal value of x 
h = increment of x used in integration
F = user-supplied function that returns the 
array F(x,y) = {y'[0],y'[1],...,y'[n-1]}.
'''

def integrate(F,x,y,xStop,h):
    
    def run_kut4(F,x,y,h):
        K0 = h*F(x,y)
        K1 = h*F(x + h/2.0, y + K0/2.0)
        K2 = h*F(x + h/2.0, y + K1/2.0)
        K3 = h*F(x + h, y + K2)
        return (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
    
    X =[]
    Y =[]
    X.append(x)
    Y.append(y)
    while x < xStop:
        h = min(h,xStop - x)
        y = y + run_kut4(F,x,y,h)
        x = x + h
        X.append(x)
        Y.append(y)
    return np.array(X),np.array(Y)

It works fine for other differential equations.

In this case the function F is defined as :

And the main code is :

def F(x,y):
    F = np.zeros(2)
    F[0] = y[1]
    F[1] = (2*(1-x)/x**3)*y[0]**(-1/2)
    return F

x = 1.0
xStop = 20
y = np.array([0,0])
h = 0.2
X,Y = integrate(F,x,y,xStop,h)
plt.plot(X,Y)
plt.grid()
plt.show()

Unfortunately, I got this error :

<ipython-input-8-8216949e6888>:4: RuntimeWarning: divide by zero encountered in power
  F[1] = (2*(1-x)/x**3)*y[0]**(-1/2)
<ipython-input-8-8216949e6888>:4: RuntimeWarning: divide by zero encountered in double_scalars
  F[1] = (2*(1-x)/x**3)*y[0]**(-1/2)

It's related to the fact that the initial value of the function is 0 but I don't know how to get rid of it in order to simplify the problem again...

Could someone help me to find an other alternative ?

Thank you for your help,

Answer 1

you y is [0,0] and in y[0]**(-1/2) there is division operation with 0 in the denominator which is giving ZeroDivision warning and invalid value encountered in double_scalars is due to expression y[0]**(-1/2) changed to NaN. however, those are warnings and F is returning value array([ 0., nan]). you need to replace y[0]**(-1/2) as negative powers of zero are undefined or you can use an extremely small value near zero if it suits your need. maybe your equation is not continuous at (1,0).

Answer 2

You can test approximate solutions close to the initial point that are powers of (1-x) or power series in this difference. In the simplest case you get y(x)=(1-x)^2 for x <= 1. To get solutions for x>1 you need to take the other sign in the square root. In total this gives a far-field behavior of x(t)=1+c*exp(-t).

Now you can follow two strategies,

• integrate the reduced equation from some initial point x=1-h with y(1-h)=h^2, y'(1-h)=-2h, along with a time integration dt/dx=y(x)^(-1/2) where t(1-h) is arbitrary, or
• integrate the original equation from some time t=T with conditions following the derivatives of x(t)=1+c*exp(T-t), that is, x(T)=1+c with arbitrary small c, x'(t)=-c, x''(T)=c. In theory they should be the same, practically with fixed-step RK4 in both cases there will be differences.

---

Make a cut across all approaches. Close to the asymptote x(t)=1 the solution is monotonous and thus time can be expressed in terms of x-1. This means that the derivative can (probably) be expressed as power series in x-1

x'  (t) = c_1 * (x-1) + c_2 * (x-1)^2 + ...
x'' (t) = c_1 * x'(t) + 2c_2 * (x-1)*x' + ...
        = (c_1 + 2c_2*(x-1)+...)*(c_1+c_2*(x-1)+..)*(x-1)
        = c_1^2*(x-1)+3c_1c_2*(x-1)^2 + ...
x'''(t) = (c_1^2+6c_1c_2*(x-1)+...)*(c_1+c_2*(x-1)+..)*(x-1)
        = c_1^3*(x-1) + 7c_1^2c_2*(x-1)^2 + ...

and 

(1-x)/x^3 = -(x-1)*(1+(x-1))^(-3)=-(x-1)+3*(x-1)^2 + ...

so equating coefficients

c_1^3=-1 ==> c_1 = -1
7c_1^2c_2 = 3 ==> c_2 = 3/7

Given x(T) close enough to 1, the other initial values have to be 
x'(T)=-(x(T)-1) + 3/7*(x(T)-1)^2
x''(T)=x(T)-1 -9/7*(x(T)-1)^2

Then the far-field approximation due to these first two terms is the solution to

x'(t) = -(x-1) + 3/7 * (x-1)^2

Substitute u(t) = (x-1)^(-1) - 3/7

u'(t) = u(t)

(x(t)-1)^(-1) - 3/7 = ((x(T)-1)^(-1) - 3/7) * exp(t-T)

x(t) = 1 + (x(T)-1)*exp(T-t) / ( 1 - 3/7*(x(T)-1)*(1-exp(T-t)) )