布拉休斯方程如下：

假設二階導數(shù)作為一階導數(shù)的斜率赚爵，做1階插值法迭代計算一階導數(shù)的值棉胀，叫做一階

這是一個非線性常微分方程，下面我們利用四階龍格庫塔方法來求解該方程冀膝，相當于找了3個值去插入斜率乘以微元量唁奢。

我們引入新的變量：

則四階布拉休斯方程可以等價的表示如下：

具有五階精度
利用四階龍格庫塔有：

對于（3）式中的第一式利用龍哥庫塔方法，有：

同理窝剖，對于(3)中的第二式：

對于(3)中的第三式：

可執(zhí)行代碼

import numpy

from matplotlib import pyplot

def GetEtaFun( eta, f1, f2, alpha ):
    eeta = []
    ff1  = []
    ff2  = []
    vv   = []
    n = len( eta )
    a13 = alpha **( 1.0/ 3.0 )
    a23 = alpha **( 2.0/ 3.0 )
    oa13 = 1.0 / a13
    for i in range( n ):
        my_eta = eta[ i ] * oa13
        my_f1 = a13 * f1[ i ]
        my_f2 = a23 * f2[ i ]
        
        v = 0.5 * ( my_eta * my_f2 - my_f1 )

        ff1.append( my_f1 )
        ff2.append( my_f2 )
        eeta.append( my_eta )
        vv.append( v )
    return eeta, ff1, ff2, vv

def RungeKutta():
    deta = 0.01
    eta  = 0

    etaN = 10.0
    n = int( etaN/deta )
    print("n=",n)
    
    f1 = 0.0
    f2 = 0.0
    f3 = 1.0
    ff1 =[f1]
    ff2 =[f2]
    ff3 =[f3]
    eeta=[eta]
    
    half = 0.5
    h = deta;
    y1 = [f1, f2, f3 ]
    y2 = [0, 0, 0]
    y3 = [0, 0, 0]
    y4 = [0, 0, 0]
    k1 = [0, 0, 0]
    k2 = [0, 0, 0]
    k3 = [0, 0, 0]
    k4 = [0, 0, 0]
    
    for i in range( n ):
        k1[0] = y1[1]
        k1[1] = y1[2]
        k1[2] = - half * y1[0] * y1[2]
       
        y2[0] = y1[0] + half * h * k1[0]
        y2[1] = y1[1] + half * h * k1[1]
        y2[2] = y1[2] + half * h * k1[2]
        
        k2[0] = y2[1]
        k2[1] = y2[2]
        k2[2] = - half * y2[0] * y2[2]
        
        y3[0] = y1[0] + half * h * k2[0]
        y3[1] = y1[1] + half * h * k2[1]
        y3[2] = y1[2] + half * h * k2[2]
        
        k3[0] = y3[1]
        k3[1] = y3[2]
        k3[2] = - half * y3[0] * y3[2]
        
        y4[0] = y1[0] + h * k3[0]
        y4[1] = y1[1] + h * k3[1]
        y4[2] = y1[2] + h * k3[2]
        
        k4[0] = y4[1]
        k4[1] = y4[2]
        k4[2] = - half * y4[0] * y4[2]        
        
        coef = 1.0 / 6.0 * h
       
        y1[0] = y1[0] + coef * ( k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0] )
        y1[1] = y1[1] + coef * ( k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1] )
        y1[2] = y1[2] + coef * ( k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2] )
        
        
        eta += h
        
        print( eta, y1[0], y1[1], y1[2] )
        ff1.append(y1[0])
        ff2.append(y1[1])
        ff3.append(y1[2])
        eeta.append(eta)
    alpha = y1[1]**(-3/2.0)
    print( "alpha=", alpha, "nstep=", n )
    
    etaList, f1List, f2List, vList = GetEtaFun( eeta, ff1, ff2, alpha )
    
    #pyplot.plot( eeta, ff1 )
    #pyplot.plot( eeta, ff2 )
    #pyplot.plot( eeta, ff3 )
    pyplot.figure(1)
    pyplot.plot( etaList, f2List )
    pyplot.figure(2)
    pyplot.plot( etaList, vList )
    pyplot.figure(3)
    pyplot.plot( f2List, etaList )
    pyplot.figure(4)
    pyplot.plot( vList, etaList )
    pyplot.show()


def main():
    RungeKutta()
    
if __name__ == "__main__":
    main()

執(zhí)行發(fā)現(xiàn)結(jié)果接近0.332