How to solve an ODE in Excel/VBA using Runge-Kutta without creating separate columns for each function evaluation.

Technorati Tags: Numerical Methods,Programming,Excel

Consider an Excel sheet to solve the compound pendulum equation numerically

The equation itself does not matter. What is important is that the angular acceleration depends on the angle and the angular speed (as well as possibly time).

The values m, L, g etc. are constants and can be declared as named ranges (for example cell C2 contains the value for m and is declared with the name "m").

I setup a row on Excel with the initial conditions, and the calculation above in a formula as the last cell in the row. For each calculation row, the position and velocity is calculated using a VBA script and called from the Excel sheet

Row	Time (C)	θ (D)	θp (E)	θpp (F)	Step (G)
9	0	0	0	=L*m*COS(D9)*(L*E9^2*SIN(D9)+2*g)/(4*I+m*L^2*COS(D9)^2)	0.02
10	=C9+G9	=RK4STEPX(G9,C9,D9,E9,F9)	=RK4STEPV(G9,C9,D9,E9,F9)	=L*m*COS(D10)*(L*E10^2*SIN(D10)+2*g)/(4*I+m*L^2*COS(D10)^2)	0.02
11	=C10+G10	=RK4STEPX(G10,C10,D10,E10,F10)	=RK4STEPV(G10,C10,D10,E10,F10)	=L*m*COS(D11)*(L*E11^2*SIN(D11)+2*g)/(4*I+m*L^2*COS(D11)^2)	0.02

The trick now is to pass into RK4STEPX and RK4STEPV the formula for θpp instead of the value. To do so the function declaration is

Public Function Rk4StepX(ByVal h As Double, ByVal t As Double, ByVal q As Double, ByVal qp As Double, ByRef qpp As Range) As Double
   ...
End Function

If you notice that qpp is passed by reference as a range the formula can be extracted by qpp.FormulaR1C1. In fact the magic actually happens in a different function tasked with evaluating the acceleration based on given position and velocity (and time).

Public Function CalcAcc(ByVal t As Double, ByVal q As Double, ByVal qp As Double, ByRef qpp As Range) As Double
    Dim f As String
    f = Mid(qpp.FormulaR1C1, 2)
    f = Replace(f, "RC[-3]", CStr(t))
    f = Replace(f, "RC[-2]", CStr(q))
    f = Replace(f, "RC[-1]", CStr(qp))
    
    CalcAcc = qpp.Worksheet.Evaluate(f)
End Function

To explain the code, the cell range containing the acceleration formula is passed in qpp. The formula is extracted in f as a string, with the leading "=" removed. The the value in the 3rd column to the left "RC[-3]" is replaced by the given value of time, the 3nd column to the left "RC[-2]" with position, and the column to the left "RC[-1]" with velocity. Lastly, the function string is evaluated on the worksheet where any named ranges exist. The expression will result in a number if the cell it corresponds resolves to a number in the worksheet.

Now that the acceleration function is defined in VBA, implementing the Runge-Kutta scheme is trivial, except this is a 2nd order ODE and the typical scheme for 1st order ODE is adapted by the author.

So here is the code that advances the position and velocity in its entirety.

Public Function Rk4StepX(ByVal h As Double, ByVal t As Double, ByVal q As Double, ByVal qp As Double, ByRef qpp As Range) As Double

    Dim K0 As Double, K1 As Double, K2 As Double, K3 As Double

    Dim Q0 As Double, Q1 As Double, Q2 As Double, Q3 As Double

    Dim h2 As Double, h6 As Double, a As Double

    h2 = h / 2: h6 = h / 6

    a = qpp.Value2

    Q0 = qp + h2 * a

    K0 = CalcAcc(t + h2, q + h2 * qp, Q0, qpp)

    Q1 = qp + h2 * K0

    K1 = CalcAcc(t + h2, q + h2 * Q0, Q1, qpp)

    Q2 = qp + h2 * K1

    K2 = CalcAcc(t + h, q + h * Q1, Q2, qpp)

    Rk4StepX = q + h * (qp + h6 * (a + K0 + K1))

End Function

Public Function Rk4StepV(ByVal h As Double, ByVal t As Double, ByVal q As Double, ByVal qp As Double, ByRef qpp As Range) As Double

    Dim K0 As Double, K1 As Double, K2 As Double, K3 As Double

    Dim Q0 As Double, Q1 As Double, Q2 As Double, Q3 As Double

    Dim h2 As Double, h6 As Double, a As Double

    h2 = h / 2: h6 = h / 6

    a = qpp.Value2

    Q0 = qp + h2 * a

    K0 = CalcAcc(t + h2, q + h2 * qp, Q0, qpp)

    Q1 = qp + h2 * K0

    K1 = CalcAcc(t + h2, q + h2 * Q0, Q1, qpp)

    Q2 = qp + h2 * K1

    K2 = CalcAcc(t + h, q + h * Q1, Q2, qpp)

    Rk4StepV = qp + h6 * (a + 2 * K0 + 2 * K1 + K2)

End Function

The results are graphed (in Excel) and shown below

Methods in Rigid Body Mechanics–Part 1

OVERVIEW

In this series of posts I will demonstrate different techniques for tackling problems with one or more rigid bodies in motion. Consider the inverted pendulum problem on a cart as described in this wikipedia page.

An approximation that can be done at this stage is that each rigid body acts like a point mass, and only the forces on each center of gravity are considered. Later we can see how to use all the degrees of freedom (including rotations) to see the effects of rotational inertias and what the reactions might be at the wheels above.

Method 1 - Cartesian Coordinates

KINEMATICS

In this example there are two bodies in motion (the cart and the ball) and most simplistic approach is to consider each rigid body moving separately. The center of gravity of the cart (Point A) has two degrees of freedom $x_A$ and $y_A$ when considered free of all contacts and constraints. Similarly the center of gravity of the ball (Point B) has also two degrees of freedom $x_B$ and $y_B$. In vector form the above coordinates are $$\vec{r}_A = \begin{pmatrix} x_A \\ y_A \\ 0 \end{pmatrix}$$ $$\vec{r}_B = \begin{pmatrix} x_B \\ y_B \\ 0 \end{pmatrix}$$

The kinematics (motion description) of such a system is trivial as the velocities and accelerations are just the derivatives of the above quantities. In vector form these are $$ \vec{v}_A = \begin{pmatrix} \dot{x}_A \\ \dot{y}_A \\ 0 \end{pmatrix}$$ $$ \vec{v}_B = \begin{pmatrix} \dot{x}_B \\ \dot{y}_B \\ 0 \end{pmatrix}$$ and $$ \vec{a}_A = \begin{pmatrix} \ddot{x}_A \\ \ddot{y}_A \\ 0 \end{pmatrix}$$ $$ \vec{a}_B = \begin{pmatrix} \ddot{x}_B \\ \ddot{y}_B \\ 0 \end{pmatrix}$$

There are a total of 4 unknown degrees of freedom in the system corresponding to the 4 acceleration components of $\vec{a}_A$ and $\vec{a}_B$.

In addition, there are two constraint equations applicable here, one keeping the cart on the ground and the other keeping the ball a distance $l$ from the cart. This makes the system with two effective degrees of freedom. To include the constraint equations in the system we are going to use what is called the Pfaffian form in terms of the velocity components.

Since the cart does not move in the vertical direction this can be written as $ \hat{j}\cdot\vec{v}_A = 0 $ or in expanded form $$ \begin{pmatrix}0\\1\\0\end{pmatrix}\cdot\begin{pmatrix} \dot{x}_A \\ \dot{y}_A \\ 0 \end{pmatrix} = \dot{y}_A = 0$$ Also the ball does not move towards the cart along the rod with direction vector $\hat{n}=\begin{pmatrix}\text{-}\sin\theta & \cos\theta & 0 \end{pmatrix}^\top$ so the constraint is $\hat{n}\cdot\vec{v}_B = 0$ or in expanded form $$\begin{pmatrix}\text{-}\sin\theta \\ \cos\theta \\ 0 \end{pmatrix} \cdot \begin{pmatrix} \dot{x}_B \\ \dot{y}_B \\ 0 \end{pmatrix} = \dot{y}_B \cos\theta – \dot{x}_B \sin\theta = 0 $$

This is kind of a tricky part, because $\theta$ depends on our degrees of freedom with $$\tan\theta = \frac{x_A-x_B}{y_B-y_A}$$ In addition, to solve for the motion we need to differentiate all these expressions in order to bring them in terms of accelerations.

To summarize, this technique quickly balloons into unmanageable complexity and thus it is not used at all in practice. The closest I have seen is the computer game particle systems because all interactions between bodies are modeled with compliant forces and not Pfaffian constraints.  Next time I will simplify this technique by only considering the reduced degrees of freedom $x$ and $\theta$ and derive everything else from those.