rkqc.txt

用VB实现解常微分方程组 包括定步长四阶龙格-库塔法、自适应变步长的龙格-库塔法、改进的中点法、外推法等

Sub RKQC(Y(), DYDX(), N, X, HTRY, EPS, YSCAL(), HDID, HNEXT)
    ONE = 1#
    SAFETY = 0.9
    ERRCON = 0.0006
    FCOR = 0.066667
    Dim YTEMP(10), YSAV(10), DYSAV(10)
    PGROW = -0.2
    PSHRNK = -0.25
    XSAV = X
    For I = 1 To N
        YSAV(I) = Y(I)
        DYSAV(I) = DYDX(I)
    Next I
    H = HTRY
    Do
        HH = 0.5 * H
        Call RK4(YSAV(), DYSAV(), N, XSAV, HH, YTEMP())
        X = XSAV + HH
        Call DERIVS(X, YTEMP(), DYDX())
        Call RK4(YTEMP(), DYDX(), N, X, HH, Y())
        X = XSAV + H
        If X = XSAV Then
            Print " Stepsize not significant in RKQC."
            Exit Sub
        End If
        Call RK4(YSAV(), DYSAV(), N, XSAV, H, YTEMP())
        ERRMAX = 0#
        For I = 1 To N
            YTEMP(I) = Y(I) - YTEMP(I)
            If ERRMAX < Abs(YTEMP(I) / YSCAL(I)) Then
                ERRMAX = Abs(YTEMP(I) / YSCAL(I))
            End If
        Next I
        ERRMAX = ERRMAX / EPS
        If ERRMAX > ONE Then
            H = SAFETY * H * (ERRMAX ^ PSHRNK)
            FLAG = 1
        Else
            HDID = H
            If ERRMAX > ERRCON Then
                HNEXT = SAFETY * H * (ERRMAX ^ PGROW)
            Else
                HNEXT = 4# * H
            End If
            FLAG = 0
        End If
    Loop While FLAG = 1
    For I = 1 To N
        Y(I) = Y(I) + YTEMP(I) * FCOR
    Next I
    Erase DYSAV, YSAV, YTEMP
End Sub