bsstep.txt

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

Sub BSSTEP(Y(), DYDX(), NV, X1, HTRY, EPS, YSCAL(), HDID, HNEXT)
    IMAX = 11
    NUSE = 7
    ONE = 1#
    SHRINK = 0.95
    GROW = 1.2
    Dim YERR(10), YSAV(10), DYSAV(10), YSEQ(10), NSEQ(11)
    NSEQ(1) = 2#:   NSEQ(2) = 4#:  NSEQ(3) = 6#:  NSEQ(4) = 8#: NSEQ(5) = 12#
    NSEQ(6) = 16#:  NSEQ(7) = 24#: NSEQ(8) = 32#: NSEQ(9) = 48#
    NSEQ(10) = 64#: NSEQ(11) = 96#
    H = HTRY
    XSAV = X1
    For I = 1 To NV
        YSAV(I) = Y(I)
        DYSAV(I) = DYDX(I)
    Next I
    Do
        For I = 1 To IMAX
            Call MMID(YSAV(), DYSAV(), NV, XSAV, H, NSEQ(I), YSEQ())
            XEST = (H / NSEQ(I)) ^ 2
            Call RZEXTR(I, XEST, YSEQ(), Y(), YERR(), NV, NUSE)
            If I > 3 Then
                ERRMAX = 0#
                For J = 1 To NV
                    If Abs(YERR(J) / YSCAL(J)) > ERRMAX Then
                        ERRMAX = Abs(YERR(J) / YSCAL(J))
                    End If
                Next J
                ERRMAX = ERRMAX / EPS
                If ERRMAX < ONE Then
                    X1 = X1 + H
                    HDID = H
                    If I = NUSE Then
                        HNEXT = H * SHRINK
                    ElseIf I = NUSE - 1 Then
                        HNEXT = H * GROW
                    Else
                        HNEXT = (H * NSEQ(NUSE - 1)) / NSEQ(I)
                    End If
                    Erase NSEQ, YSEQ, DYSAV, YSAV, YERR
                    Exit Sub
                End If
            End If
        Next I
        H = 0.25 * H / 2 ^ ((IMAX - NUSE) / 2)
    Loop Until X1 + H = X1
    Print " Step size underflow"
End Sub