/*Для тестирования C#-программы принят пример, приведенный в книге Шуп. Т. 
Решение инженерных задач на ЭВМ: Практическое ру-ководство: 
Пер. с англ. – М.: Мир, 1982. – С.79 - 85. */
using System; // C#-программа
using System.IO;
/* Метод Рунге-Кутта решения системы ОДУ – 
Method Runge-Kutta decisions of the SCDE (two equations)*/
// Decision of the system of the common differential equations

namespace ConsAppl_SCDE_2
{
    class DecisionSCDE
    {
        public DecisionSCDE() { }
        private static void Main()
        {
            int n = 2;                                   // система двух ОДУ                                           
            double[] p = { 0.0, 1.0, 0.02, 0.0001 };
            double[] w = { 10.0, 0.0 };             // начальные значения искомых;

            DecisionSCDE metodRK = new DecisionSCDE();

            metodRK.SCDE(p, w, n);
        }
        private void SCDE(double[] p, double[] w, int n)
        /*Здесь реализовано решение системы ОДУ порядка n (в примере n = 2) 
          методом Рунге-Кутта с модификацией Фельберга */
        {
            double[] f = new double[n];
            double[] y = new double[n];
            double[] u = new double[n];
            double[] a = new double[n];
            double[] b = new double[n];
            double[] c = new double[n];
            double[] d = new double[n];
            double[] e = new double[n];

            double x, xk, h, e1, e2, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c13, q;
            int j, nu, klu, key;

            /* x - начало интервала интегрирования; xk - конец интервала; 
          h - начальный шаг интегрирования; e1 – точность решения системы ОДУ*/
            klu = 1; x = p[0]; xk = p[1]; h = p[2]; e1 = p[3];

            c1 = 1.0 / 9.0; c2 = 1.0 / 64.0; c3 = 1.0 / 2.4; c4 = 1.0 / 15.0; c5 = 1.0 / 144.0; c6 = 1.0 / 6.0;
            c7 = 1.0 / 45.0; c8 = 1.0 / 75.0; c9 = 1.0 / 150.0; c10 = 1.0 / 12.0; c11 = 1.0 / 432.0; c13 = 1.0 / 27.0;

            for (j = 0; j < n; j++) y[j] = w[j];
            OUTP(y, p, x, n);           //1-я печать
            loop1: key = 0;
            FCT(y, f, x);          //1. Расчет производных
            nu = 0;
            for (j = 0; j < n; j++)
            { a[j] = f[j] * h; y[j] = w[j] + 2 * a[j] * c1; }
            x += 2 * h * c1;
            FCT(y, f, x);          //2. Расчет производных

            for (j = 0; j < n; j++)
            { b[j] = h * f[j]; y[j] = w[j] + a[j] * c10 + b[j] * 0.25; }
            x += h * c1;
            FCT(y, f, x);          //3. Расчет производных
            for (j = 0; j < n; j++)
            { c[j] = f[j] * h; y[j] = w[j] + (34.5 * a[j] - 121.5 * b[j] + 135 * c[j]) * c2; }
            x += h * c3;
            FCT(y, f, x);          //4. Расчет производных
            for (j = 0; j < n; j++)
            {
                d[j] = f[j] * h; q = w[j] - 17 * a[j] * c10 + 6.75 * b[j];
                y[j] = q - 5.4 * c[j] + 16 * d[j] * c4;
            };
            x += h * 0.25;
            FCT(y, f, x);           //5. Расчет производных
            for (j = 0; j < n; j++)
            {
                e[j] = f[j] * h; q = w[j] + 65 * a[j] * c11 - 0.3125 * b[j];
                y[j] = q + 0.8125 * c[j] + 4 * d[j] * c13 + 5 * e[j] * c5;
            }
            x -= h * c6;
            FCT(y, f, x);           //6. Расчет производных
            for (j = 0; j < n; j++)
            {
                u[j] = h * f[j];
                y[j] = w[j] + a[j] * c1 + 0.45 * c[j] + 16 * d[j] * c7 + e[j] * c10;
                q = a[j] * c9 - 0.03 * c[j] + 16 * d[j] * c8 + e[j] * 0.05;
                e2 = Math.Abs(q - 0.24 * u[j]);
                if (e2 > e1) key = 1;
                else if (e2 < e1 * 0.05) nu += 1;
            }
            x += h * c6;
            if (key == 1)
            { x -= h; for (j = 0; j < n; j++) y[j] = w[j]; h *= 0.5; goto loop1; }

            if (nu == n) h += h;
            if (h > p[2]) h = p[2];
            OUTP(y, p, x, n);             //Печать
            for (j = 0; j < n; j++) w[j] = y[j];
            if (klu == 2) goto loop2;
            //калибровка последнего шага на интервале интегрирования
            if (Math.Abs(x - xk) <= h) { h = Math.Abs(x - xk); klu = 2; }
            goto loop1;
            loop2: return;
        }

        private void OUTP(double[] y, double[] p, double t, int n)
        {// печать результатов решения системы двух ОДУ (2.28)

            FileStream file = new FileStream("result.txt", FileMode.Append);
            StreamWriter write = new StreamWriter(file);
            if (t == 0)
            {
                write.WriteLine();
                write.WriteLine("Решение системы {0}-x ОДУ методом Рунге-Кутта.", n);
                write.WriteLine();
                write.WriteLine("Шаг интегрирования определяется автоматически с целью");
                write.WriteLine("достижения заданной точности.");
                write.WriteLine();
                write.WriteLine("Начало диапазона интегрирования                  tn = {0} c.", p[0]);
                write.WriteLine("Окончание диапазона интегрирования           tk = {0} c.", p[1]);
                write.WriteLine("Начальное значение шага интегрирования    dt = {0} c.", p[2]);
                write.WriteLine("Верхняя граница погрешностей                              = {0} . ", p[3]);
                write.WriteLine();
                write.WriteLine("    Результаты решения системы ОДУ: ");
                write.WriteLine(" --------------------------------------------------------");
                write.WriteLine("Время t, c \tПуть y, см \tСкорость dy/dt, см/с");
                write.WriteLine(" --------------------------------------------------------");
            }
            write.WriteLine("{0:#.######}  \t{1:#.####}  \t{2:#.####}", t, y[0], y[1]);
            if (t == 1)
            {
                write.WriteLine(" --------------------------------------------------------");
                write.WriteLine("    Программа работу закончила ");
            }
            write.Flush();
            write.Close();
            file.Close();
        }

        private double[] FCT(double[] y, double[] f, double t)
        { // вычисление значений производных: правые части  системы двух ОДУ (2.28) 
            double m_1, c, xk_1, xks;
            m_1 = 100.0; c = 0.15; xk_1 = 2.0; xks = 0.2;

            f[0] = y[1];
            f[1] = m_1 * (-c * f[0] - xk_1 * y[0] - xks * Math.Pow(y[0], 3));
            return f;
        }
    }
}
