Численные методы в задачах теплопереноса. Метод чехарда

namespace mod3._2
{
    public partial class Form1 : Form
    {
        static private double F(double x) => -(Math.Pow(x, 2) / 100) + x - 24; // Функция исходной параболы
 
        static private int N = 100;
        static private double maxX = 100;
 
        static private double hx = maxX / N;
        static private double ht = 0.5;
        static private double r = ht / Math.Pow(hx, 2);
 
        //  коэффициенты, одинаковые на каждом шаге;
        static private double a = r;
        static private double b = -1 - 2 * r;
        static private double c = r;
 
        //  коэффициенты, вычисляемые на каждом шаге
        static private double[] f = new double[N];
        static private double[] alpha = new double[N];
        static private double[] beta = new double[N];
 
        /// <summary>
        /// 
        /// Начальное рапределение температуры
        /// 
        /// </summary>
        /// <returns></returns>
        static private double[] Start()
        {
            double[] res = new double[N];
            for (int i = 0; i < N; i++)
                if (i <= 40 || i >= 60) res[i] = 0;
                else res[i] = F(i);
            return res;
        }
 
        /// <summary>
        /// 
        /// Явная схема
        /// Стр. 77
        /// 
        /// </summary>
        /// <param name="T"></param>
        /// <returns></returns>
        static private double[] Met2(double T)
        {
            double[] mas = Start();
 
            double t = 0;
            double[] present = new double[N];
            for (int i = 1; i < N - 1; i++)
                present[i] = mas[i];
 
            while (t < T)
            {
                for (int i = 1; i < N - 1; i++)
                    mas[i] = r * present[i + 1] + (1 - 2 * r) * present[i] + r * present[i - 1];
                for (int i = 1; i < N - 1; i++)
                    present[i] = mas[i];
                t += ht;
            }
 
            return mas;
        }
 
        /// <summary>
        /// 
        /// Неявная схема
        /// Стр. 78
        /// 
        /// </summary>
        /// <param name="T"></param>
        /// <returns></returns>
        static private double[] Met3(double T)
        {
            double[] mas = Start();
 
            double t = 0;
            while (t < T)
            {
                for (int i = 0; i < N; i++)
                    f[i] = -mas[i];
 
                alpha[0] = -c / b;
                beta[0] = f[0] / b;
 
                for (int i = 1; i < N - 1; i++)
                {
                    alpha[i] = -c / (a * alpha[i - 1] + b);
                    beta[i] = (f[i] - a * beta[i - 1]) / (a * alpha[i - 1] + b);
                }
 
                mas[N - 1] = (f[N - 1] - a * beta[N - 2]) / (a * alpha[N - 2] + b);
 
                for (int i = N - 2; i >= 0; i--)
                    mas[i] = alpha[i] * mas[i + 1] + beta[i];
 
                t += ht;
            }
 
            return mas;
        }
 
 
 
        public Form1()
        {
            Init();
        }
 
        /// <summary>
        /// 
        /// Отрисовка графиков
        /// 
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void Form_Load(object sender, EventArgs e)
        {
            double[] mas = new double[N];
 
            mas = Start();
            for (int i = 0; i < N; i++)
            {
                series[0].Points.AddXY(i * hx, mas[i]);
                series[0].Points[i].BorderWidth = 3;
            }
            mas = Met2(10);
            for (int i = 0; i < N; i++)
            {
                series[1].Points.AddXY(i * hx, mas[i]);
                series[1].Points[i].BorderWidth = 3;
            }
            mas = Met2(25);
            for (int i = 0; i < N; i++)
            {
                series[2].Points.AddXY(i * hx, mas[i]);
                series[2].Points[i].BorderWidth = 3;
            }
            mas = Met2(50);
            for (int i = 0; i < N; i++)
            {
                series[3].Points.AddXY(i * hx, mas[i]);
                series[3].Points[i].BorderWidth = 3;
            }
            mas = Met2(75);
            for (int i = 0; i < N; i++)
            {
                series[4].Points.AddXY(i * hx, mas[i]);
                series[4].Points[i].BorderWidth = 3;
            }
            mas = Met2(100);
            for (int i = 0; i < N; i++)
            {
                series[5].Points.AddXY(i * hx, mas[i]);
                series[5].Points[i].BorderWidth = 3;
            }
 
 
 
            mas = Start();
            for (int i = 0; i < N; i++)
            {
                series[6].Points.AddXY(i * hx, mas[i]);
                series[6].Points[i].BorderWidth = 3;
            }
            mas = Met3(10);
            for (int i = 0; i < N; i++)
            {
                series[7].Points.AddXY(i * hx, mas[i]);
                series[7].Points[i].BorderWidth = 3;
            }
            mas = Met3(25);
            for (int i = 0; i < N; i++)
            {
                series[8].Points.AddXY(i * hx, mas[i]);
                series[8].Points[i].BorderWidth = 3;
            }
            mas = Met3(50);
            for (int i = 0; i < N; i++)
            {
                series[9].Points.AddXY(i * hx, mas[i]);
                series[9].Points[i].BorderWidth = 3;
            }
            mas = Met3(75);
            for (int i = 0; i < N; i++)
            {
                series[10].Points.AddXY(i * hx, mas[i]);
                series[10].Points[i].BorderWidth = 3;
            }
            mas = Met3(100);
            for (int i = 0; i < N; i++)
            {
                series[11].Points.AddXY(i * hx, mas[i]);
                series[11].Points[i].BorderWidth = 3;
            }
        }
    }
}