Python ODE system 4th order with equation system

Question

I'm coding a python script for the solution for the following ode system. I use scipy solve_bvp library.

Here is my implementation:

from scipy.integrate import solve_bvp
import numpy as np


def ode(r, y):
    Re=100
    return np.vstack((
#       y[0]    #  f        
        y[1],   # [f]' = [f']
        2/r**y[0]-2/3*Re/r**2*y[3]*y[0],   #[f']'=f''= ...
#       y[3]    #  g        
        y[4],   # [g]'=g'
        y[5],   # [g']'=g''
        y[6],   # [g'']'=g'''
        12/r**2*y[4]-24/r**3*y[4]+2*Re/(5*r**3)*(
            5*y[0]*(r*y[1]-2*y[0])+
            2*r*y[3]*y[6]+
            (r*y[4]-4*y[3])*y[5]
            -18/r*y[4]*y[5]
            +48/r**2*y[3]**2
            )
    ))


def bc(ya, yb):
    return np.array([
                     ya[0],     yb[0]-1,
                     ya[1],     yb[1],
                     ya[3],     yb[3], 
                     ya[4],     yb[4], 
                     ya[5],     yb[5],
                     ya[6],     yb[6]
                     ])

x_flow = np.linspace(0.7, 1, 10)
y_flow = np.ones((6, x_flow.shape[0]))

res_flow = solve_bvp(ode, bc, x_flow, y_flow, verbose=2)

Currently the solver return the error:

  y[6],   # [g'']'=g'''

IndexError: index 6 is out of bounds for axis 0 with size 6

Using the debugger, the error is thrown on line 647 of _bvp.py while wrapping the custom ode with np.asarray().

Answer 1

You have 6 variables: [f0, f0', g1, g1', g1'', g1'''] - this should be reflected in your y variable, your derivative vector and your boundary conditions, noting that the corresponding indices go from [0] to [5].

This looks like a similarity solution for something - sorry, I haven't worked out which - so although the following code runs it may or may not produce what you want. Note that f0 is nearly linear at this Reynolds number and radius range, so isn't plotted below.

It's also useful to start from a guess that satisfies the boundary conditions at least.

from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
import numpy as np

Re = 100.0                   # Reynolds number
eta = 0.7                    # inner radius

def ode( r, y ):             # y is [ f0, f0', g1, g1', g1'', g1''' ]
    [ f0, f0p, g1, g1p, g1pp, g1ppp ] = y
    RHSf = ( 2 * f0 / r ** 2 ) * ( 1 - Re * g1p / 3 ),
    RHSg = ( 12 * g1pp / r ** 2 - 24 * g1p / r ** 3 
         + ( 2 * Re / 5 / r ** 3 ) * (   5 * f0 * ( r * f0p - 2 * f0 )
                                       + 2 * r * g1 * g1ppp
                                       + ( r * g1p - 4 * g1 ) * g1pp - 18 * g1 * g1p / r + 48 * g1 ** 2 / r ** 2 ) )
    return np.vstack( ( f0p, RHSf, g1p, g1pp, g1ppp, RHSg ) )


def bc( ya, yb ):
    return np.array( [ ya[0], ya[2], ya[3],      yb[0] - 1.0, yb[2], yb[3] ] )


r = np.linspace( eta, 1.0, 100 )
y = np.zeros( ( 6, len( r ) ) )
y[0] = np.linspace( 0.0, 1.0, len( r ) )        # initial guess to fit boundary conditions

result = solve_bvp( ode, bc, r, y, verbose=2 )
rplot = np.linspace( eta, 1.0, 100 )
f = result.sol( rplot )[0]
g = result.sol( rplot )[2]

#plt.plot( r, f, label="f" )
plt.plot( r, g, label="g" )
plt.legend()
plt.xlabel( "r" )
plt.show()

Output (for g1):

Answer 2

The issue is with the indexes. You are creating a y vector with only 6 elements (so the index can go from 0 to 5) but then you are calling y[6] so the 7th element. However just creating a y_flow with 7 elements like

y_flow = np.ones((7, x_flow.shape[0]))

does not work because you have only 6 values returned by your ODE, so either you have error in the indexing of y in the ode function or you are missing an equation. I added a dummy equation uncommenting the first y[0] in the ODE and got an error in the boundary condition because you are passing 12 values instead of 7.

I think you need to reformulate your system a bit better before putting it into code in order to be consistent with the indexing.

My suggestion is to write accurately every equation in the same form as you would write in the code and then code it and the same goes for the boundary conditions.