odeint.txt

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

Sub ODEINT(YSTART(), NVAR, X1, X2, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
    MAXSTP = 10000
    TWO = 2#
    ZERO = 0#
    TINY = 1E-30
    Dim YSCAL(10), Y(10), DYDX(10)
    X = X1
    H = Abs(H1) * Sgn(X2 - X1)
    NOK = 0
    NBAD = 0
    KOUNT = 0
    For I = 1 To NVAR
        Y(I) = YSTART(I)
    Next I
    If KMAX > 0 Then XSAV = X - DXSAV * TWO
    For NSTP = 1 To MAXSTP
        Call DERIVS(X, Y(), DYDX())
        For I = 1 To NVAR
            YSCAL(I) = Abs(Y(I)) + Abs(H * DYDX(I)) + TINY
        Next I
        If KMAX > 0 Then
            If Abs(X - XSAV) > Abs(DXSAV) Then
                If KOUNT < KMAX - 1 Then
                    KOUNT = KOUNT + 1
                    XP(KOUNT) = X
                    For I = 1 To NVAR
                        YP(I, KOUNT) = Y(I)
                    Next I
                    XSAV = X
                End If
            End If
        End If
        If (X + H - X2) * (X + H - X1) > ZERO Then H = X2 - X
        Call RKQC(Y(), DYDX(), NVAR, X, H, EPS, YSCAL(), HDID, HNEXT)
        If HDID = H Then
            NOK = NOK + 1
        Else
            NBAD = NBAD + 1
        End If
        If (X - X2) * (X2 - X1) >= ZERO Then
            For I = 1 To NVAR
                YSTART(I) = Y(I)
            Next I
            If KMAX <> 0 Then
                KOUNT = KOUNT + 1
                XP(KOUNT) = X
                For I = 1 To NVAR
                    YP(I, KOUNT) = Y(I)
                Next I
            End If
            Erase DYDX, Y, YSCAL
            Exit Sub
        End If
        If Abs(HNEXT) < HMIN Then
            Print " Stepsize smaller than minimum."
            Exit Sub
        End If
        H = HNEXT
    Next NSTP
    Print " Too many steps."
End Sub