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,

Is this from your previous `f'''(x)=(1-f)/f^3`, reducing the order by setting `f'(x)=u(f(x))` to get `f''(x)=u'(f(x))*u(f(x))` and `f'''(x)=u''(f(x))*u(f(x))^2+u'(f(x))^2*u(f(x))`? Then the equation becomes `u''(f)=-u'(f)^2/u(f)+(1-f)/f^3/u(f)^2`. Is that on your path or did you try a different substitution? — Lutz Lehmann, Dec 14 '20 at 15:40
It is Carl Runge and Martin Wilhelm Kutta, two very "kaiserlich-deutsche" names. — Lutz Lehmann, Dec 14 '20 at 15:43
Thank you Lutz Lehmann for the reference : yes it is from my previous post. In fact, I used a slightly different substitution by setting : `u(f(x)) = f ' (x)^2` so we have `u '' (f) = 2*f '''(x) * u**(-1/2)`. Due to the fact that the boundary conditions of the previous problem were at infinity : f(inf) = 1, f ' (inf) = 0 and f ' ' (inf) = 0. I think that the new problem seems to be non continuous at 1.... — Wiss, Dec 14 '20 at 16:44
This parametrization is only valid for strictly monotonous solutions, or restrictions of solutions to such segments. The reference should be the constant solution `f(x)=1` that is not covered by the parametrization. Close to the constant solution you have in the linearization of the difference to it terms `exp(-x)` that are indeed falling towards this constant, but also terms `exp(x/2)*cos(sqrt(3)*x/2)` and the sine variant, which are not monotonous and increasing in amplitude. This alone makes the numerical solution rather unstable, or the problem ill-posed. — Lutz Lehmann, Dec 14 '20 at 17:09
Thanks a lot for your explanations : I have tried the substitution that you have proposed earlier setting `f'(x)=u(f(x))`. I got a solution that doesn't seem coherent with what I have expected since I have implemented initial conditions at values different from zero but very close to... It seems to be the best solution version for the moment.... — Wiss, Dec 14 '20 at 20:11

hack3r-0m · Answer 1 · 2020-12-14T14:50:31.203

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).

edited Dec 14 '20 at 14:50

answered Dec 14 '20 at 13:44

hack3r-0m

700
6
20

thank you for your answer : indeed it was what I've suspected but I'm sorry but I didn't get your advice when you said "check parentheses..." ? could you specify this point please ? thank you again – Wiss Dec 14 '20 at 13:56
parentheses seem to be correct as compared to the pic, I have updated my answer – hack3r-0m Dec 14 '20 at 14:37
Thank you for your precisions : I tried to take a small value near zero and calculations seem to be done until x = 1.2... after that an other error raises : `RuntimeWarning: invalid value encountered in double_scalars F[1] = (2*(1-x)/x**3)*y[0]**(-1/2)` I don't understand how to manipulate the boundary conditions again in order to simplify this equation if it is not continuous... thank you again – Wiss Dec 14 '20 at 14:57

Lutz Lehmann · Answer 2 · 2020-12-15T11:13:25.990

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)) )

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

2 Answers2