DotNumerics Implicit Runge Kutta

Linear Algebra

Linear Equations.

Optimization

Simplex.

Differential Equations

Implicit Runge Kutta method.

The OdeImplicitRungeKutta5 class solves an initial-value problem for stiff ordinary differential equations using the implicit Runge Kutta method of order 5.

In this example the OdeImplicitRungeKutta5 class is used to solve the ODEs for a kinetic model:

              a1       a2
         S0 <---> S1 <---> S2
              b1       b2
The differential equations are:
            y0' = -y0 * a1 + y1 * b1;
            y1' = -y1 * (b1 + a2) + y0 * a1 + y2 * b2;
            y2' = -y2 * b2 + y1 * a2;
where
        a1 = 0.01* Exp(0.7 * voltage / 25)
        a2 = 0.01* Exp(0.4 * voltage / 25)

C# Code

using System;

using System.Collections.Generic;

using System.Windows.Forms;

using DotNumerics.ODE;

namespace DotNumericsDemo.DifferentialEquations

    public partial class ImplicitRungeKutta5 : Form

        private OdeImplicitRungeKutta5 implicitRK5 = new OdeImplicitRungeKutta5();

        private double[] t = new double[1001];

        private double[] v = new double[1001];

        private double[] S0 = new double[1001];

        private double[] S1 = new double[1001];

        private double[] S2 = new double[1001];

        private int solIndex = 0;

        private double voltage = 0;

        private double voltage1 = 50;

        private double voltage2 = -150;

        private double voltage3 = 50;

        private void Solve()

            OdeFunction fun = new OdeFunction(ODEs);

            OdeSolution sol = new OdeSolution(OdeSol);

            this.implicitRK5.InitializeODEs(fun, 3);

            solIndex = 0;

            double[] y0 = new double[3];

            //Solve for voltage1

            y0[0] = 1;

            y0[1] = 0;

            y0[2] = 0;

            this.voltage = this.voltage1;

            this.implicitRK5.Solve(y0, 0, 1, 300, sol);

            //Solve for voltage2

            y0[0] = this.S0[300];

            y0[1] = this.S1[300];

            y0[2] = this.S2[300];

            this.voltage = this.voltage2;

            this.implicitRK5.Solve(y0, 301, 1, 700, sol);

            //Solve for voltage3

            y0[0] = this.S0[700];

            y0[1] = this.S1[700];

            y0[2] = this.S2[700];

            this.voltage = this.voltage3;

            this.implicitRK5.Solve(y0, 701, 1, 1000, sol);

        //      a1       a2

        // S0 <---> S1 <---> S2

        //      b1       b2

        //Differential equations:

        //            y'0 = -y0 * a1 + y1 * b1;

        //            y'1 = -y1 * (b1 + a2) + y0 * a1 + y2 * b2;

        //            y'2 = -y2 * b2 + y1 * a2;

        //where

        //        a1 = 0.01* Exp(0.7 * voltage / 25)

        //        a2 = 0.01* Exp(0.4 * voltage / 25)

        double[] yprime = new double[3];

        private double[] ODEs(double t, double[] y)

            double a1 = 0.01 * Math.Exp(0.7 * this.voltage / 25);

            double a2 = 0.01 * Math.Exp(0.4 * this.voltage / 25);

            double b1 = 0.001 * Math.Exp(-0.7 * this.voltage / 25);

            double b2 = 0.001 * Math.Exp(-0.4 * this.voltage / 25);

            yprime[0] = -y[0] * a1 + y[1] * b1;

            yprime[1] = -y[1] * (b1 + a2) + y[0] * a1 + y[2] * b2;

            yprime[2] = -y[2] * b2 + y[1] * a2;

            return this.yprime;

        //Save the solution

        private void OdeSol(double t, double[] y)

            this.t[solIndex] = t;

            this.v[solIndex] = voltage;

            this.S0[solIndex] = y[0];

            this.S1[solIndex] = y[1];

            this.S2[solIndex] = y[2];

            solIndex++;