pythonodesolver

Python ODE system 4th order with equation system

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

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

Solution

  • 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): enter image description here