' these are the initial values used in the primary
' example of section 8.6 of Boyce & DiPrima'a
' "Elementary Differential Equations and Boundary Value Problems"

' we have an ODE of the form
'     y' = 1 - x + 4 * y
' we want to solve for values of y over the interval from
' x = 0 to x = 1.  When x = 0, y = 1.
' We are going to break this up into 100 subintervals.

' first, we get a RungeKutta4 instance:

set rk4 = New RungeKutta4

' it doesn't matter in what order we set the properties;
' we just need to get them all set before we solve.

' first we'll tell it what equation we are solving:
rk4.ODE = "1 - x + 4 * y"

' next, we enter the initial conditions (x=0, y=1).  Note that
' our "x" is "t" in the RungeKutta class:
rk4.t0 = 0: rk4.y0 = 1

' ..and we will be solving for y as x goes to 1:

rk4.tmax = 1

' Finally we want to test over 100 intervals, so:
rk4.intervals = 100

' solving the equation returns an array of solution points in s
s = rk4.solve

' we can join each one for easier display
for i = 0 to UBound(s)
        s(i) = Join(s(i), vbTab)
next

' then do another join to allow single-call output of the table
s = Join(s, vbLf)
WScript.Echo s




Class RungeKutta4

        ' a VBScript Implementation of a 4th-order Runge-Kutta solution to an
        ' ordinary differential equation
        ' based on reference from Section 8.6 of Boyce & DiPrima's 4th Edition
        ' "Elementary Differential Equations and Boundary Value Problems"
        '        
        ' Basic Runge-Kutta implementation
        ' provide the following values for the calculation:
        ' + initial values (t0, y0)
        ' + final t value
        ' + number of steps between t0 and t
        ' + ODE - the first order differential to be solved.
        '
        ' the calculate method will return an array with an element
        ' for each step from t0 to t; each element will consist of
        ' the calculated t and y for the step.
        '
        ' CAVEATS ON USE:
        ' All numerical methods have limitations.  It is best to ensure continuity
        ' of the range over the solution domain when using them; also, since
        ' numerical methods will always exhibit rounding errors, at a certain
        ' point using smaller step sizes to decrease discretization errors will
        ' cause even greater roundoff discussion of the issues involved with this.


        dim m_t0, m_y0, m_h, m_tmax, m_intervals, m_ymax


        Property Let t0(val)
                m_t0 = val
        End Property


        Property Let y0(val)
                m_y0 = val
        End Property


        Property get ymax
                ' found in calculate
                ' this is the value of y at the end of the interval of interest;
                ' useful if you don't need intermediary results on the interval.
                ymax = m_ymax
        end property


        Property Let ODE(equation)
                dim fn
                fn = Array("function RHS(t,y)", "RHS = CDbl( " & equation & " )",  _
                        "End function")
                ExecuteGlobal Join(fn, vbLf)
        end property


        Property Let tmax(val)
                m_tmax = val
        End Property


        property let intervals(num)
                m_intervals = num
                m_h = (m_tmax - m_t0)/num
        end property


        function solve
                dim rk(), i, kk(4)
                redim rk(m_intervals)
                t = m_t0
                y = m_y0
                rk(0) = Array(m_t0, m_y0)
                for i = 1 to m_intervals
                        kk(1) = RHS(t, y)
                        kk(2) = RHS(t + .5*m_h, y + .5*m_h*kk(1))
                        kk(3) = RHS(t + .5*m_h, y + .5*m_h*kk(2))
                        kk(4) = RHS(t + m_h, y + m_h*kk(3))
                        y = y + m_h/6 * ( kk(1) + 2*kk(2) + 2*kk(3) + kk(4) )
                        t = m_t0 + (m_tmax - m_t0)/m_intervals*i
                        rk(i) = Array(t,y)
                next
                m_ymax = y
                solve = rk
        end function


End Class