Решение уравнений методом Рунге-Кутта - C#

 class Program
    {
        public delegate double Function(double x, double y);
 
        #region Ordinary Differential Equations - Methods
 
        public static double ODE_Euler(Function f, double x0, double y0, double h, double x)
        {
            double xnew, ynew, result = double.NaN;
            if (x <= x0)
                result = y0;
            else if (x > x0)
            {
                do
                {
                    if (h > x - x0) h = x - x0;
                    ynew = y0 + f(x0, y0) * h;
                    xnew = x0 + h;
                    x0 = xnew;
                    y0 = ynew;
                } while (x0 < x);
                result = ynew;
            }
            return result;
        }
 
        public static double ODE_RungeKutta2(Function f, double x0, double y0, double h, double x)
        {
            double xnew, ynew, k1, k2, result = double.NaN;
            if (x == x0)
                result = y0;
            else if (x > x0)
            {
                do
                {
                    if (h > x - x0) h = x - x0;
                    k1 = h * f(x0, y0);
                    k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
                    ynew = y0 + k2;
                    xnew = x0 + h;
                    x0 = xnew;
                    y0 = ynew;
                } while (x0 < x);
                result = ynew;
            }
            return result;
        }
 
        public static double ODE_RungeKutta4(Function f, double x0, double y0, double h, double x)
        {
            double xnew, ynew, k1, k2, k3, k4, result = double.NaN;
            if (x == x0)
                result = y0;
            else if (x > x0)
            {
                do
                {
                    if (h > x - x0) h = x - x0;
                    k1 = h * f(x0, y0);
                    k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
                    k3 = h * f(x0 + 0.5 * h, y0 + 0.5 * k2);
                    k4 = h * f(x0 + h, y0 + k3);
                    ynew = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
                    xnew = x0 + h;
                    x0 = xnew;
                    y0 = ynew;                                        
                }while (x0 < x);
                result = ynew;
            }
            return result;
        }
 
        public static double ODE_RungeKuttaFehlberg(Function f, double x0, double y0, double x, double h, double tolerance)
        {
            double xnew, ynew, hnew, k1, k2, k3, k4, k5, k6;
            double hmin = 0.0001;
            double hmax = 0.5;
            if (h > hmax) h = hmax;
            if (h < hmin) h = hmin;
 
            while (x0 < x)
            {
                k1 = h * f(x0, y0);
                k2 = h * f(x0 + 0.25 * h, y0 + 0.25 * k1);
                k3 = h * f(x0 + 3 * h / 8, y0 + 3 * k1 / 32 + 9 * k2 / 32);
                k4 = h * f(x0 + 12 * h / 13, y0 + 1932 * k1 / 2197 - 7200 * k2 / 2197 + 7296 * k3 / 2197);
                k5 = h * f(x0 + h, y0 + 439 * k1 / 216 - 8 * k2 + 3680 * k3 / 513 - 845 * k4 / 4104);
                k6 = h * f(x0 + 0.5 * h, y0 - 8 * k1 / 27 + 2 * k2 - 3544 * k3 / 2565 + 1859 * k4 / 4104 - 11 * k5 / 40);
                double error = Math.Abs(k1 / 360 - 128 * k3 / 4275 - 2197 * k4 / 75240 + k5 / 50 + 2 * k6 / 55) / h;
                double s = Math.Pow(0.5 * tolerance / error, 0.25);
                if (error < tolerance)
                {
                    ynew = y0 + 25 * k1 / 216 + 1408 * k3 / 2565 + 2197 * k4 / 4104 - 0.2 * k5;
                    xnew = x0 + h;
                    x0 = xnew;
                    y0 = ynew; 
                }
                if (s < 0.1) s = 0.1;
                if (s > 4)   s = 4;
                hnew = h*s;
                h = hnew;
                if (h > hmax) h = hmax;
                if (h < hmin) h = hmin;
                if (h > x - x0) h = x - x0;
            } return y0;
        }
 
        #endregion
 
        #region Ordinary Differential Equations - Test functions
 
        static double f(double x, double y)
        {
            return y*Math.Cos(x);
        }
 
        static double dx(double x, double y, double z)
        { return 10.0 * (y - x); }
 
        static double dy(double x, double y, double z)
        { return x * (28.0 - z) - y; }
 
        static double dz(double x, double y, double z)
        { return x * y - 8.0 * z / 3.0; }
 
        static void TestEuler()
        {
            double h = 0.001;
            double x0 = 0.0;
            double y0 = 1.0;
            Console.WriteLine("\n Results from Euler's method with h = {0}\n", h);
            double result = y0;
            for (int i = 0; i < 11; i++)
            {
                double x = 0.1 * i;
                result = ODE_Euler(f, x0, result, h, x);
                double exact = Math.Exp(Math.Sin(x));
                if (i % 5 == 0)
                Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
                x0 = x;
            }
        }
        
        static void TestRungeKutta2()
        {
            double h = 0.001;
            double x0 = 0.0;
            double y0 = 1.0;
            Console.WriteLine("\n Results from the 2nd-order Runge-Kutta method with h = {0}\n", h);
            double result = y0;
            for (int i = 0; i < 11; i++)
            {
                double x = 0.1 * i;
                result = ODE_RungeKutta2(f, x0, result, h, x);
                double exact = Math.Exp(Math.Sin(x));
                if (i % 5 == 0)
                Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
                x0 = x;
            }
        }
 
        static void TestRungeKutta4()
        {
            double h = 0.001;
            double x0 = 0.0;
            double y0 = 1.0;
            Console.WriteLine("\n Results from the 4th-order Runge-Kutta method with h = {0}\n", h);
            double result = y0;
            for (int i = 0; i < 11; i++)
            {
                double x = 0.1 * i;
                result = ODE_RungeKutta4(f, x0, result, h, x);
                double exact = Math.Exp(Math.Sin(x));
                if (i % 5 == 0)
                Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
                x0 = x;
            }
        }
 
        static void TestRungeKuttaFehlberg()
        {
            double h = 0.2;
            double x0 = 0.0;
            double y0 = 1.0;
            Console.WriteLine("\n Results from the fourth-order Runge-Kutta-Fehlberg method with h = {0}\n", h);
            double result = y0;
            for (int i = 0; i < 11; i++)
            {
                double x = 0.1 * i;
                result = ODE_RungeKuttaFehlberg(f, x0, result, x, h, 1e-8);
                double exact = Math.Exp(Math.Sin(x));
                if (i%5==0)
                Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
                x0 = x;
            }
        }
 
        static void TestEulerCoupledODE()
        {
            double x = 0.0, y = 5.0, z = 10.0;  //initial conditions
            double dt = 0.001;                  //step size
 
            Console.WriteLine("\n Results from a coupled ODE using Euler's method\n");
 
            for (int i = 0; i < 100; i++)
            {
                double xnew = x + dx(x, y, z) * dt;
                double ynew = y + dy(x, y, z) * dt;
                double znew = z + dz(x, y, z) * dt;
                x = xnew;
                y = ynew;
                z = znew;
                Console.WriteLine(" x = {0:n1}, y = {1:e12}, z = {2:e12}", x, y, z);
            }
        }
        #endregion
        
        static void Main(string[] args)
        {
            TestEuler();
            TestRungeKutta2();
            TestRungeKutta4();
            TestRungeKuttaFehlberg();
            TestEulerCoupledODE();
            Console.ReadLine();
        }
    }