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

How to solve an elastic spring system with C#

How to write a simple FEA solver in C#. There is only one type of element, a linear spring. In the example below there are 5 nodes and 4 elements. They form a tree structure one branch connects nodes 1 through 4, and the second branch connects nodes 3 and 5. There is spring pre-loading between elements 2 and 3, and an external force applied on node 5. Also the stiffness varies between elements. Nodes 1 and 4 are immobile (fixed to the ground). The constraints and loadings have been highlighted in bold below.

The program result is:

NODE LIST
#         POS      DEFL     FORCE       FIX
1       0.000   0.00000     1.304     FIXED
2       1.000  -0.13043     0.000      FREE
3       2.000  -0.04348     0.000      FREE
4       3.000   0.00000     8.696     FIXED
5       3.000  -0.14348   -10.000      FREE

ELEMENT LIST
#   NODES STIFFNESS   PRELOAD     FORCE
1     1→2    10.000     0.000     1.304
2     2→3   100.000    10.000     1.304
3     3→4   200.000     0.000    -8.696
4     3→5   100.000     0.000    10.000

The full code listing is below:
/// <summary>
/// Nodes can either have prescribed force 
/// or displacement.
/// </summary>
public enum Apply
{
    Force,
    Displacement
}

public class Node
{
    /// <summary>
    /// Define position x, 
    /// deformation u, 
    /// and applied force W
    /// </summary>
    public double x, ux, W;
    /// <summary>
    /// Define if displacement or force 
    /// is applied to node
    /// </summary>
    public Apply fix;

    /// <summary>
    /// Defines a free node at a position.
    /// </summary>
    /// <param name="x"></param>
    public Node(double x) : this(x, Apply.Force, 0) { }
    /// <summary>
    /// Defines a node at a position, 
    /// with prescribed force or displacement.
    /// </summary>
    /// <param name="x">The node position</param>
    /// <param name="fix">constraint type</param>
    /// <param name="value">constraint value</param>
    public Node(double x, Apply fix, double value)
    {
        this.x=x;
        this.fix=fix;
        if(fix==Apply.Force)
        {
            this.W=value;
            this.ux=0;
        }
        else
        {
            this.ux=value;
            this.W=0;
        }
    }
}

public class Element
{
    /// <summary>
    /// Connects node i with node j
    /// </summary>
    public int i, j;
    /// <summary>
    /// Stiffness k, preload F0
    /// and element force f
    /// </summary>
    public double k, F0, f;

    /// <summary>
    /// Define new element between two nodes
    /// </summary>
    /// <param name="i">From node index</param>
    /// <param name="j">To node index</param>
    /// <param name="k">Element stiffness</param>
    public Element(int i, int j, double k)
        : this(i, j, k, 0) { }
    /// <summary>
    /// Define new element between two nodes
    /// </summary>
    /// <param name="i">From node index</param>
    /// <param name="j">To node index</param>
    /// <param name="k">Element stiffness</param>
    /// <param name="F0">Element pre-load</param>
    public Element(int i, int j, double k, double F0)
    {
        this.i=i;
        this.j=j;
        this.k=k;
        this.F0=F0;
        this.f=0;
    }
}

/// <summary>
/// Create a linear system of springs
/// </summary>
public class WorldSystem
{
    public List<Node> nodes;
    public List<Element> elements;
    Matrix P, A, R1, R2;
    SymmetricBandMatrix K;

    public WorldSystem()
    {
        this.nodes=new List<Node>();
        this.elements=new List<Element>();
    }

    #region Structure
    /// <summary>
    /// Add a free node
    /// </summary>
    /// <param name="x">The node position</param>
    public void AddFreeNode(double x)
    {
        nodes.Add(new Node(x));
    }
    /// <summary>
    /// Add a free node with force
    /// </summary>
    /// <param name="x">The node position</param>
    /// <param name="W">The applied force</param>
    public void AddFreeNode(double x, double W)
    {
        nodes.Add(new Node(x, Apply.Force, W));
    }
    /// <summary>
    /// Add a fixed node
    /// </summary>
    /// <param name="x">The node position</param>
    public void AddFixedNode(double x)
    {
        nodes.Add(new Node(x, Apply.Displacement, 0));
    }
    /// <summary>
    /// Add a fixed node with displacement
    /// </summary>
    /// <param name="x">The node position</param>
    /// <param name="ux">The applied displacement</param>
    public void AddFixedNode(double x, double ux)
    {
        nodes.Add(new Node(x, Apply.Displacement, ux));
    }
    /// <summary>
    /// Add a spring between two nodes
    /// </summary>
    /// <param name="i">From node index</param>
    /// <param name="j">To node index</param>
    /// <param name="k">Spring stiffness</param>
    /// <param name="F0">Spring preload</param>
    public void AddElement(int i, int j,
        double k, double F0)
    {
        elements.Add(new Element(i, j, k, F0));
    }
    /// <summary>
    /// Add a spring between two nodes
    /// </summary>
    /// <param name="i">From node index</param>
    /// <param name="j">To node index</param>
    /// <param name="k">Spring stiffness</param>
    public void AddElement(int i, int j, double k)
    {
        elements.Add(new Element(i, j, k));
    }
    /// <summary>
    /// Add a spring between two nodes
    /// </summary>
    /// <param name="A">From node</param>
    /// <param name="B">To node</param>
    /// <param name="k">Spring stiffness</param>
    public void AddElement(Node A, Node B, double k)
    {
        int i=nodes.IndexOf(A);
        int j=nodes.IndexOf(B);
        AddElement(i, j, k);
    }
    /// <summary>
    /// Add a spring between two nodes
    /// </summary>
    /// <param name="A">From node</param>
    /// <param name="B">To node</param>
    /// <param name="k">Spring stiffness</param>
    /// <param name="F0">Spring preload</param>
    public void AddElement(Node A, Node B,
        double k, double F0)
    {
        int i=nodes.IndexOf(A);
        int j=nodes.IndexOf(B);
        AddElement(i, j, k, F0);
    }
    #endregion

    /// <summary>
    /// Create an example system with 
    /// 5 nodes and 4 elements
    /// </summary>
    public static WorldSystem CreateTestCase()
    {
        //                  v==3==(D)
        // (A)==1==(B)==2==(C)
        //                  ^==4==(E)
        //
        // Fixed nodes (A), (D)
        // Force applied on node (E)

        WorldSystem sys=new WorldSystem();

        // Fixed Node (A)
        sys.AddFixedNode(0.0);
        // Free Node (B)
        sys.AddFreeNode(1.0);
        // Free Node (C)
        sys.AddFreeNode(2.0);
        // Fixed (D)
        sys.AddFixedNode(3.0);
        // Free Node (E) with force applied
        sys.AddFreeNode(3.0, -10.0);

        // Element (A)===(B)
        sys.AddElement(0, 1, 10);
        // Element (B)===(C)
        sys.AddElement(1, 2, 100, 10);
        // Element (C)===(D)
        sys.AddElement(2, 3, 200);
        // Element (C)===(E)
        sys.AddElement(2, 4, 100);

        return sys;
    }


    #region Solution
    public void SolveSystem(ref Vector b, ref Vector x)
    {
        // Split vectors x & b into two components
        //
        // x = R1*x1 + R2*x2
        // where x1=known(H), x2=unknown(N+M-H)
        //
        // b = R1*b1 + R2*b2
        // where b1=unknown(H), b2=known(N+M-H)

        // Filter known nodal displacements
        var x1=R1.Transpose()*x;
        // Filter for known element preloads
        // and nodal forces
        var b2=R2.Transpose()*b;

        var lse=new LinearEquations();

        // reduce the system by substituting for known
        // displacements for forces
        var A22=R2.Transpose()*A*R2;
        var b22=R2.Transpose()*(R2*b2-A*R1*x1);
        // solve for unknown element forces
        // and nodal displacements
        var x2=lse.Solve(A22, b22);
        x=(R1*x1+R2*x2).GetColumnVector(0);
        // back substitute for unknown nodal forces
        var b1=R1.Transpose()*A*x;
        b=(R1*b1+R2*b2).GetColumnVector(0);
    }

    public void ResetSystem()
    {
        Vector b, x;
        BuildSolution(out b, out x);
        // Filter out unknown quantities
        x=(R1*R1.Transpose()*x).GetColumnVector(0);
        b=(R2*R2.Transpose()*b).GetColumnVector(0);
        ReportSolution(b, x);
    }

    void BuildSystem()
    {
        int N=nodes.Count;
        int M=elements.Count;

        // Connection Matrix P:{#ELEM, #NODES}
        this.P=new Matrix(M, N);
        // Stiffness Matrix K:diag{#ELEM}
        this.K=new SymmetricBandMatrix(M, 0);

        this.A=new Matrix(M+N, M+N);
        for(int k=0; k<M; k++)
        {
            P[k, elements[k].i]=1;
            P[k, elements[k].j]=-1;
            K[k, k]=elements[k].k;
        }
        var PT=P.Transpose();
        var KP=K*P;
        // Element Forces:      f = F - K*P*u
        // Nodal Forces:        W - PT*f = 0

        // As a system of equations A*x=b
        //
        // | [1] -[K]*[P] | | [f] |   | [F] |
        // |              | |     | = |     |
        // | [P]T    0    | | [u] |   | [W] |
        //
        // At this stage some of [u] may be known 
        // and some of [W] may be unknown
        for(int k=0; k<M; k++)
        {
            A[k, k]=1;
        }
        int H=0; // Count fixed nodes
        for(int i=0; i<N; i++)
        {
            if(nodes[i].fix==Apply.Displacement) H++;
            for(int k=0; k<M; k++)
            {
                A[k, M+i]=-KP[k, i];
                A[M+i, k]=PT[i, k];
            }
        }
        // Build the known and unknown projection matrices

        // Example with 3 nodes & 2 elements. 
        // 2 nodes have prescribed displacements
        //
        // | f1 |   | 0 0 |          | 1 0 0 |        
        // | f2 |   | 0 0 | | u1 |   | 0 1 0 | | f1 | 
        // | u1 | = | 1 0 | |    | + | 0 0 0 | | f2 | 
        // | u2 |   | 0 0 | | u3 |   | 0 0 1 | | u2 | 
        // | u3 |   | 0 1 |          | 0 0 0 |        
        //
        //  x = R1*x1 + R2*x2

        this.R1=new Matrix(N+M, H);
        this.R2=new Matrix(N+M, N+M-H);
        for(int k=0; k<M; k++)
        {
            // all of internal forces are unknown
            R2[k, k]=1;
        }
        for(int i=0, j=0, l=0; i<N; i++)
        {
            // set projection matrix coefficients
            if(nodes[i].fix==Apply.Displacement)
            {
                R1[M+i, j++]=1;
            }
            else
            {
                R2[M+i, M+l++]=1;
            }
        }
    }
    void BuildSolution(out Vector b, out Vector x)
    {
        int N=nodes.Count;
        int M=elements.Count;

        b=new Vector(M+N);
        x=new Vector(M+N);
        for(int k=0; k<M; k++)
        {
            // get element forces into solution 
            // Vectors b, x
            b[k]=elements[k].F0;
            x[k]=elements[k].f;
        }
        for(int i=0; i<N; i++)
        {
            // get nodal forces and displacements
            // into Vectors b, x
            b[M+i]=nodes[i].W;
            x[M+i]=nodes[i].ux;
        }
    }
    void ReportSolution(Vector b, Vector x)
    {
        int N=nodes.Count;
        int M=elements.Count;

        for(int k=0; k<M; k++)
        {
            // set element forces from solution
            // Vectors b, x
            elements[k].f=x[k];
            elements[k].F0=b[k];
        }
        for(int i=0; i<N; i++)
        {
            // set nodal forces and displacements
            // from Vectors b, x
            nodes[i].ux=x[M+i];
            nodes[i].W=b[M+i];
        }
    }
    #endregion

    #region Report
    public override string ToString()
    {
        StringBuilder sb=new StringBuilder();
        sb.AppendLine("NODE LIST");
        sb.AppendFormat(
            "{0,-3} {1,9} {2,9} {3,9} {4,9}",
            "#",
            "POS",
            "DEFL",
            "FORCE",
            "FIX");
        sb.AppendLine();
        for(int i=0; i<nodes.Count; i++)
        {
            sb.AppendFormat(
                "{0,-3} {1,9:F3} {2,9:F5} {3,9:F3} {4,9}",
                i+1,
                nodes[i].x,
                nodes[i].ux,
                nodes[i].W,
                nodes[i].fix==Apply.Displacement?"FIXED":"FREE");
            sb.AppendLine();
        }

        sb.AppendLine("ELEMENT LIST");
        sb.AppendFormat(
            "{0,-3} {1,5} {2,9} {3,9} {4,9}",
            "#",
            "NODES",
            "STIFFNESS",
            "PRELOAD",
            "FORCE");
        sb.AppendLine();
        for(int k=0; k<elements.Count; k++)
        {
            var nodes=(elements[k].i+1).ToString()+"→"+(elements[k].j+1).ToString();
            sb.AppendFormat(
                "{0,-3} {1,5} {2,9:F3} {3,9:F3} {4,9:F3}",
                k+1,
                nodes,
                elements[k].k,
                elements[k].F0,
                elements[k].f);
            sb.AppendLine();

        }
        return sb.ToString();
    }

    #endregion

    static void Main(string[] args)
    {
        var lss=CreateTestCase();
        Vector b, x;
        // Build system strucutre
        lss.BuildSystem();
        // Move knowns from system to b, x
        lss.BuildSolution(out b, out x);
        // Solve for the unknowns
        lss.SolveSystem(ref b, ref x);
        // Move unknowns into system
        lss.ReportSolution(b, x);

        // Report results
        Console.WriteLine(lss.ToString());
        Console.WriteLine("Press ENTER to quit.");
        Console.ReadLine();
    }
}

Technorati Tags: Programming,FEM,C#,Engineering