fork download
copy 
  1. using System;
  2.  
  3. public class Test
  4. {
  5. 	public static void Main()
  6. 	{
  7. 		// your code goes here
  8. 	}
  9. 	public interface LinearAlgebra
  10. {
  11.     private const double EPS = 1e-12;
  12.     private const double TINY = 1e-20;
  13.  
  14.     private const double DEFAULT_EPS = 1e-6;
  15.     private const double SECOND_DIFF_EPS = 1e-4;
  16.  
  17.     private static bool IsEqual(double a, double b) => Math.Abs(a - b) < EPS;
  18.  
  19.     // Глубокое копирование jagged-массива
  20.     private static double[][] DeepCopy(double[][] original)
  21.     {
  22.         if (original == null) return null;
  23.         double[][] copy = new double[original.Length][];
  24.         for (int i = 0; i < original.Length; i++)
  25.         {
  26.             if (original[i] != null)
  27.             {
  28.                 copy[i] = new double[original[i].Length];
  29.                 Array.Copy(original[i], copy[i], original[i].Length);
  30.             }
  31.         }
  32.         return copy;
  33.     }
  34.  
  35.     // Копирование вектора
  36.     private static double[] DeepCopy(double[] original)
  37.     {
  38.         if (original == null) return null;
  39.         double[] copy = new double[original.Length];
  40.         Array.Copy(original, copy, original.Length);
  41.         return copy;
  42.     }
  43.  
  44.     /// <summary>
  45.     /// Решает систему A * x = b методом LU-разложения.
  46.     /// Исходные матрица A и вектор b НЕ изменяются.
  47.     /// </summary>
  48.     /// <param name="a">Квадратная матрица коэффициентов (jagged array)</param>
  49.     /// <param name="b">Вектор правой части</param>
  50.     /// <returns>Вектор решения x</returns>
  51.     public static double[] SolveLinearSystem(double[][] a, double[] b)
  52.     {
  53.         if (a == null || b == null)
  54.             throw new ArgumentNullException();
  55.         int n = a.Length;
  56.         if (n == 0 || b.Length != n)
  57.             throw new ArgumentException("Matrix must be square and vector size must match n.");
  58.         for (int i = 0; i < n; i++)
  59.         {
  60.             if (a[i] == null || a[i].Length != n)
  61.                 throw new ArgumentException($"Row {i} is null or has incorrect length.");
  62.         }
  63.  
  64.         // Работаем с копиями!
  65.         double[][] aCopy = DeepCopy(a);
  66.         double[] bCopy = DeepCopy(b);
  67.  
  68.         int[] indx = new int[n];
  69.         double d = 1.0;
  70.  
  71.         LU_Decomposition(ref aCopy, indx, ref d, n);
  72.         LU_Solve(aCopy, ref bCopy, indx, n);
  73.  
  74.         return bCopy; // теперь это решение x
  75.     }
  76.  
  77.     private static void LU_Decomposition(ref double[][] a, int[] indx, ref double d, int n)
  78.     {
  79.         double[] vv = new double[n];
  80.         d = 1.0;
  81.  
  82.         for (int i = 0; i < n; i++)
  83.         {
  84.             double big = 0.0;
  85.             for (int j = 0; j < n; j++)
  86.             {
  87.                 double temp = Math.Abs(a[i][j]);
  88.                 if (temp > big) big = temp;
  89.             }
  90.             if (IsEqual(big, 0.0))
  91.                 throw new InvalidOperationException("Singular matrix in LU_Decomposition.");
  92.             vv[i] = 1.0 / big;
  93.         }
  94.  
  95.         for (int j = 0; j < n; j++)
  96.         {
  97.             for (int i = 0; i < j; i++)
  98.             {
  99.                 double sum = a[i][j];
  100.                 for (int k = 0; k < i; k++)
  101.                     sum -= a[i][k] * a[k][j];
  102.                 a[i][j] = sum;
  103.             }
  104.  
  105.             double big = 0.0;
  106.             int imax = j;
  107.             for (int i = j; i < n; i++)
  108.             {
  109.                 double sum = a[i][j];
  110.                 for (int k = 0; k < j; k++)
  111.                     sum -= a[i][k] * a[k][j];
  112.                 a[i][j] = sum;
  113.  
  114.                 double dum = vv[i] * Math.Abs(sum);
  115.                 if (dum > big || IsEqual(dum, big))
  116.                 {
  117.                     big = dum;
  118.                     imax = i;
  119.                 }
  120.             }
  121.  
  122.             if (j != imax)
  123.             {
  124.                 (a[imax], a[j]) = (a[j], a[imax]);
  125.                 d = -d;
  126.                 vv[imax] = vv[j];
  127.             }
  128.  
  129.             indx[j] = imax;
  130.  
  131.             if (IsEqual(a[j][j], 0.0))
  132.                 a[j][j] = TINY;
  133.  
  134.             if (j != n - 1)
  135.             {
  136.                 double dum = 1.0 / a[j][j];
  137.                 for (int i = j + 1; i < n; i++)
  138.                     a[i][j] *= dum;
  139.             }
  140.         }
  141.     }
  142.  
  143.     private static void LU_Solve(double[][] a, ref double[] b, int[] indx, int n)
  144.     {
  145.         int ii = -1;
  146.  
  147.         for (int i = 0; i < n; i++)
  148.         {
  149.             int ip = indx[i];
  150.             double sum = b[ip];
  151.             b[ip] = b[i];
  152.  
  153.             if (ii != -1)
  154.             {
  155.                 for (int j = ii; j < i; j++)
  156.                     sum -= a[i][j] * b[j];
  157.             }
  158.             else if (!IsEqual(sum, 0.0))
  159.             {
  160.                 ii = i;
  161.             }
  162.             b[i] = sum;
  163.         }
  164.  
  165.         for (int i = n - 1; i >= 0; i--)
  166.         {
  167.             double sum = b[i];
  168.             for (int j = i + 1; j < n; j++)
  169.                 sum -= a[i][j] * b[j];
  170.  
  171.             double ld = a[i][i];
  172.             b[i] = IsEqual(ld, 0.0) ? 0.0 : sum / ld;
  173.         }
  174.     }
  175.  
  176.  
  177.     /// <summary>
  178.     /// Вычисляет обратную матрицу для заданной квадратной матрицы.
  179.     /// Исходная матрица НЕ изменяется.
  180.     /// </summary>
  181.     /// <param name="a">Исходная квадратная матрица (jagged array)</param>
  182.     /// <returns>Обратная матрица A⁻¹</returns>
  183.     public static double[][] InvertMatrix(double[][] a)
  184.     {
  185.         if (a == null)
  186.             throw new ArgumentNullException(nameof(a));
  187.         int n = a.Length;
  188.         if (n == 0)
  189.             throw new ArgumentException("Matrix cannot be empty.");
  190.         for (int i = 0; i < n; i++)
  191.         {
  192.             if (a[i] == null || a[i].Length != n)
  193.                 throw new ArgumentException($"Row {i} is null or has incorrect length.");
  194.         }
  195.  
  196.         // Создаём единичную матрицу
  197.         double[][] identity = CreateIdentityMatrix(n);
  198.         double[][] inverse = new double[n][];
  199.  
  200.         // Для каждого столбца единичной матрицы решаем A * x = e_j
  201.         for (int j = 0; j < n; j++)
  202.         {
  203.             double[] e_j = new double[n];
  204.             for (int i = 0; i < n; i++)
  205.                 e_j[i] = identity[i][j]; // j-й столбец
  206.  
  207.             // Решаем систему A * x = e_j
  208.             double[] x_j = SolveLinearSystem(a, e_j);
  209.  
  210.             // x_j — это j-й столбец обратной матрицы
  211.             // Сохраняем его как j-й столбец в inverse
  212.             if (inverse[0] == null) // инициализируем строки, если ещё не созданы
  213.             {
  214.                 for (int i = 0; i < n; i++)
  215.                     inverse[i] = new double[n];
  216.             }
  217.  
  218.             for (int i = 0; i < n; i++)
  219.                 inverse[i][j] = x_j[i];
  220.         }
  221.  
  222.         return inverse;
  223.     }
  224.  
  225.     // Вспомогательный метод: создаёт единичную матрицу n×n
  226.     private static double[][] CreateIdentityMatrix(int n)
  227.     {
  228.         double[][] I = new double[n][];
  229.         for (int i = 0; i < n; i++)
  230.         {
  231.             I[i] = new double[n];
  232.             I[i][i] = 1.0;
  233.         }
  234.         return I;
  235.     }
  236.  
  237.  
  238.  
  239.     /// <summary>
  240.     /// Умножает матрицу A (m×n) на матрицу B (n×p): C = A * B
  241.     /// </summary>
  242.     public static double[][] MultiplyMatrices(double[][] a, double[][] b)
  243.     {
  244.         if (a == null || b == null)
  245.             throw new ArgumentNullException();
  246.         if (a.Length == 0 || b.Length == 0)
  247.             throw new ArgumentException("Matrices cannot be empty.");
  248.  
  249.         int m = a.Length;           // строки A
  250.         int n = a[0].Length;        // столбцы A (и строки B)
  251.         int p = b[0].Length;        // столбцы B
  252.  
  253.         // Проверка совместимости размеров
  254.         if (b.Length != n)
  255.             throw new ArgumentException("Number of columns in A must equal number of rows in B.");
  256.  
  257.         // Проверка, что все строки A и B имеют правильную длину
  258.         for (int i = 0; i < m; i++)
  259.             if (a[i] == null || a[i].Length != n)
  260.                 throw new ArgumentException($"Row {i} of matrix A is invalid.");
  261.         for (int i = 0; i < n; i++)
  262.             if (b[i] == null || b[i].Length != p)
  263.                 throw new ArgumentException($"Row {i} of matrix B is invalid.");
  264.  
  265.         // Результат: m×p
  266.         double[][] result = new double[m][];
  267.         for (int i = 0; i < m; i++)
  268.         {
  269.             result[i] = new double[p];
  270.             for (int j = 0; j < p; j++)
  271.             {
  272.                 double sum = 0.0;
  273.                 for (int k = 0; k < n; k++)
  274.                 {
  275.                     sum += a[i][k] * b[k][j];
  276.                 }
  277.                 result[i][j] = sum;
  278.             }
  279.         }
  280.  
  281.         return result;
  282.     }
  283.  
  284.     /// <summary>
  285.     /// Умножает матрицу A (m×n) на вектор v (n): y = A * v
  286.     /// </summary>
  287.     public static double[] MultiplyMatrixByVector(double[][] a, double[] v)
  288.     {
  289.         if (a == null || v == null)
  290.             throw new ArgumentNullException();
  291.         if (a.Length == 0 || v.Length == 0)
  292.             throw new ArgumentException("Matrix or vector cannot be empty.");
  293.  
  294.         int m = a.Length;
  295.         int n = v.Length;
  296.  
  297.         // Проверка: все строки A должны иметь длину n
  298.         for (int i = 0; i < m; i++)
  299.         {
  300.             if (a[i] == null || a[i].Length != n)
  301.                 throw new ArgumentException($"Row {i} of matrix A does not match vector length.");
  302.         }
  303.  
  304.         double[] result = new double[m];
  305.         for (int i = 0; i < m; i++)
  306.         {
  307.             double sum = 0.0;
  308.             for (int k = 0; k < n; k++)
  309.             {
  310.                 sum += a[i][k] * v[k];
  311.             }
  312.             result[i] = sum;
  313.         }
  314.  
  315.         return result;
  316.     }
  317.  
  318.  
  319.     /// <summary>
  320.     /// Вычисляет определитель квадратной матрицы с использованием LU-разложения.
  321.     /// </summary>
  322.     /// <param name="a">Квадратная матрица (jagged array)</param>
  323.     /// <returns>Определитель матрицы det(A)</returns>
  324.     public static double Determinant(double[][] a)
  325.     {
  326.         if (a == null)
  327.             throw new ArgumentNullException(nameof(a));
  328.         int n = a.Length;
  329.         if (n == 0)
  330.             throw new ArgumentException("Matrix cannot be empty.");
  331.         for (int i = 0; i < n; i++)
  332.         {
  333.             if (a[i] == null || a[i].Length != n)
  334.                 throw new ArgumentException($"Row {i} is null or has incorrect length.");
  335.         }
  336.  
  337.         // Работаем с копией, чтобы не менять оригинал
  338.         double[][] aCopy = DeepCopy(a);
  339.         int[] indx = new int[n];
  340.         double d = 1.0;
  341.  
  342.         // Выполняем LU-разложение (с pivoting)
  343.         LU_Decomposition(ref aCopy, indx, ref d, n);
  344.  
  345.         // Определитель = d * произведение диагональных элементов U
  346.         double det = d;
  347.         for (int i = 0; i < n; i++)
  348.         {
  349.             det *= aCopy[i][i];
  350.         }
  351.  
  352.         return det;
  353.     }
  354.  
  355.     /// <summary>
  356.     /// Вычисляет C = A + alpha * B (поэлементно)
  357.     /// </summary>
  358.     /// <param name="a">Вектор A</param>
  359.     /// <param name="b">Вектор B</param>
  360.     /// <param name="alpha">Коэффициент</param>
  361.     /// <returns>Новый вектор C</returns>
  362.     public static double[] AddVectorsWithScale(double[] a, double[] b, double alpha)
  363.     {
  364.         if (a == null || b == null)
  365.             throw new ArgumentNullException();
  366.         if (a.Length != b.Length)
  367.             throw new ArgumentException("Vectors must have the same length.");
  368.  
  369.         double[] c = new double[a.Length];
  370.         for (int i = 0; i < a.Length; i++)
  371.         {
  372.             c[i] = a[i] + alpha * b[i];
  373.         }
  374.         return c;
  375.     }
  376.  
  377.  
  378.  
  379.     /// <summary>
  380.     /// Численно вычисляет градиент скалярной функции R: ℝⁿ → ℝ в точке x.
  381.     /// </summary>
  382.     /// <param name="R">Функция R(x)</param>
  383.     /// <param name="x">Точка, в которой вычисляется градиент</param>
  384.     /// <param name="eps">Шаг конечных разностей (по умолчанию 1e-8)</param>
  385.     /// <returns>Вектор градиента ∇R(x) размера n</returns>
  386.     public static double[] Gradient(Func<double[], double> R, double[] x, double eps = DEFAULT_EPS)
  387.     {
  388.         if (R == null) throw new ArgumentNullException(nameof(R));
  389.         if (x == null) throw new ArgumentNullException(nameof(x));
  390.  
  391.         int n = x.Length;
  392.         double[] grad = new double[n];
  393.         double[] xPerturbed = new double[n];
  394.         Array.Copy(x, xPerturbed, n);
  395.  
  396.         double f0 = R(x); // f(x)
  397.  
  398.         for (int i = 0; i < n; i++)
  399.         {
  400.             double original = x[i];
  401.             xPerturbed[i] = original + eps;
  402.             double f1 = R(xPerturbed);
  403.             grad[i] = (f1 - f0) / eps;
  404.             xPerturbed[i] = original; // восстановить
  405.         }
  406.  
  407.         return grad;
  408.     }
  409.  
  410.     /// <summary>
  411.     /// Численно вычисляет якобиан вектор-функции R: ℝⁿ → ℝᵐ в точке x.
  412.     /// Якобиан J — матрица m×n, где J[i,j] = ∂R_i/∂x_j
  413.     /// </summary>
  414.     /// <param name="R">Вектор-функция R(x) → double[m]</param>
  415.     /// <param name="x">Точка в ℝⁿ</param>
  416.     /// <param name="eps">Шаг конечных разностей</param>
  417.     /// <returns>Якобиан: J[i][j]</returns>
  418.     public static double[][] Jacobian(Func<double[], double[]> R, double[] x, double eps = DEFAULT_EPS)
  419.     {
  420.         if (R == null) throw new ArgumentNullException(nameof(R));
  421.         if (x == null) throw new ArgumentNullException(nameof(x));
  422.  
  423.         double[] y0 = R(x);
  424.         if (y0 == null) throw new InvalidOperationException("R(x) returned null.");
  425.         int m = y0.Length;
  426.         int n = x.Length;
  427.  
  428.         double[][] J = new double[m][];
  429.         for (int i = 0; i < m; i++)
  430.             J[i] = new double[n];
  431.  
  432.         double[] xPerturbed = new double[n];
  433.         Array.Copy(x, xPerturbed, n);
  434.  
  435.         for (int j = 0; j < n; j++)
  436.         {
  437.             double original = x[j];
  438.             xPerturbed[j] = original + eps;
  439.             double[] y1 = R(xPerturbed);
  440.  
  441.             for (int i = 0; i < m; i++)
  442.             {
  443.                 J[i][j] = (y1[i] - y0[i]) / eps;
  444.             }
  445.  
  446.             xPerturbed[j] = original;
  447.         }
  448.  
  449.         return J;
  450.     }
  451.  
  452.  
  453.     /// <summary>
  454.     /// Численно вычисляет гессиан скалярной функции R: ℝⁿ → ℝ в точке x.
  455.     /// Гессиан H — симметричная матрица n×n, где H[i,j] = ∂²R/∂x_i∂x_j
  456.     /// </summary>
  457.     /// <param name="R">Скалярная функция R(x)</param>
  458.     /// <param name="x">Точка в ℝⁿ</param>
  459.     /// <param name="eps">Шаг конечных разностей</param>
  460.     /// <returns>Гессиан: H[i][j]</returns>
  461.  
  462.     public static double[][] Hessian(Func<double[], double> R, double[] x, double eps = SECOND_DIFF_EPS)
  463.     {
  464.         if (R == null) throw new ArgumentNullException(nameof(R));
  465.         if (x == null) throw new ArgumentNullException(nameof(x));
  466.  
  467.         int n = x.Length;
  468.         double[][] H = new double[n][];
  469.         for (int i = 0; i < n; i++)
  470.             H[i] = new double[n];
  471.  
  472.         double[] xTemp = new double[n];
  473.         Array.Copy(x, xTemp, n);
  474.  
  475.         double f0 = R(x); // f(x)
  476.  
  477.         for (int i = 0; i < n; i++)
  478.         {
  479.             for (int j = i; j < n; j++)
  480.             {
  481.                 if (i == j)
  482.                 {
  483.                     // Вторая производная по x_i
  484.                     double orig = x[i];
  485.                     xTemp[i] = orig + eps;
  486.                     double f1 = R(xTemp);        // f(x + h e_i)
  487.                     xTemp[i] = orig + 2 * eps;
  488.                     double f2 = R(xTemp);        // f(x + 2h e_i)
  489.                     xTemp[i] = orig;             // восстановить
  490.  
  491.                     H[i][i] = (f2 - 2 * f1 + f0) / (eps * eps);
  492.                 }
  493.                 else
  494.                 {
  495.                     // Смешанная производная ∂²f/∂x_i∂x_j
  496.                     double origI = x[i], origJ = x[j];
  497.  
  498.                     // f(x + h e_i + h e_j)
  499.                     xTemp[i] = origI + eps;
  500.                     xTemp[j] = origJ + eps;
  501.                     double fPP = R(xTemp);
  502.  
  503.                     // f(x + h e_i)
  504.                     xTemp[j] = origJ;
  505.                     double fP0 = R(xTemp);
  506.  
  507.                     // f(x + h e_j)
  508.                     xTemp[i] = origI;
  509.                     xTemp[j] = origJ + eps;
  510.                     double f0P = R(xTemp);
  511.  
  512.                     // f(x)
  513.                     xTemp[i] = origI;
  514.                     xTemp[j] = origJ;
  515.                     // f0 уже известен
  516.  
  517.                     H[i][j] = (fPP - fP0 - f0P + f0) / (eps * eps);
  518.                     H[j][i] = H[i][j]; // симметрия
  519.                 }
  520.             }
  521.         }
  522.  
  523.         return H;
  524.     }
  525.  
  526.     /// <summary>
  527.     /// Возвращает полный гиперкуб, полученный путем перебора значений переменных 
  528.     /// от заданного минимального значения до заданного максимального значения.
  529.     /// По всем переменным одинаковое количество точек от максимума до минимума.        
  530.     /// </summary>
  531.     /// <param name="variablesMinMax">Двухмерный массив. который содержит значения минимального и максимального значения для каждой переменной</param>
  532.     /// <param name="pointsPerVariables">Количество перебираемых точек для каждой переменной</param>
  533.     /// <returns>Гиперкуб, заполненный всеми возможными комбинациями. Размер растет экспоненциально от количества точек</returns>
  534.     /// <exception cref="ArgumentNullException"></exception>
  535.     static double[][] bruteForceGrid(double[][] variablesMinMax, int pointsPerVariables)
  536.     {
  537.         if (variablesMinMax == null) throw new ArgumentNullException(nameof(variablesMinMax));
  538.         if (pointsPerVariables < 2) pointsPerVariables = 2;
  539.  
  540.         int Variables = (int) variablesMinMax.Length;
  541.         long M = (long) Math.Pow(pointsPerVariables, Variables);
  542.         double [][] totalCombinationsArray = new double[M][];
  543.         int[] indexes = new int[Variables];
  544.  
  545.         for (long i = 0; i < M; i++)
  546.         {
  547.             long c = i;
  548.             int varCounter = 0;
  549.  
  550.             totalCombinationsArray[i] = new double[Variables];
  551.  
  552.             for (varCounter = 0; varCounter < Variables; varCounter++)
  553.             {
  554.                 indexes[varCounter] = (int) c % pointsPerVariables;
  555.                 c /= pointsPerVariables;
  556.  
  557.                 double min = variablesMinMax[varCounter][0];
  558.                 double max = variablesMinMax[varCounter][1];
  559.                 double delta = (max - min) / (pointsPerVariables-1) * indexes[varCounter];
  560.                 double value = min + delta;
  561.                 totalCombinationsArray[i][varCounter] = value;
  562.             }   
  563.  
  564.         }
  565.  
  566.         return totalCombinationsArray;
  567.  
  568.     }
  569.  
  570.     /// <summary>
  571.     /// Находит минимум функции R(x) полным перебором по сетке (полному гиперкубу). 
  572.     /// Полный гиперкуб, получается путем перебора значений переменных
  573.     /// от заданного минимального значения до заданного максимального значения.
  574.     /// По всем переменным одинаковое количество точек от максимума до минимума.        
  575.     /// </summary>
  576.     /// <param name="variablesMinMax">Двухмерный массив. который содержит значения минимального и максимального значения для каждой переменной</param>
  577.     /// <param name="pointsPerVariables">Количество перебираемых точек для каждой переменной</param>
  578.     /// <returns>Вектор x, при котором достигает минимум функции R(x)</returns>        
  579.  
  580.     public static double[] BruteForceMinimize(Func<double[], double> R, double[][] variablesMinMax, int pointsPerVariables = 5)
  581.     {
  582.         int N = variablesMinMax.Length;
  583.         double[][] totalCombinationsArray = bruteForceGrid(variablesMinMax, pointsPerVariables);
  584.         long Variants = totalCombinationsArray.LongLength;
  585.         double[] X = new double[N];
  586.         X = totalCombinationsArray[0];
  587.         double Ri = R(X), Rmin;
  588.         Rmin = Ri;
  589.  
  590.         for (long i = 0; i < Variants; i++)
  591.         {
  592.             Ri = R(totalCombinationsArray[i]);
  593.             if (Ri < Rmin && Ri > 0)
  594.             {
  595.                 Rmin = Ri;
  596.                 X = totalCombinationsArray[i];
  597.             }
  598.         }
  599.  
  600.         return X;
  601.  
  602.     }
  603.  
  604.  
  605.  
  606.     /// <summary>
  607.     /// Минимизирует скалярную функцию R(x) методом Ньютона.
  608.     /// </summary>
  609.     /// <param name="R">Функция для минимизации: R(x) → double</param>
  610.     /// /// <param name="Validator">Функция для проверки вектора x на физичность: Validator(x) → bool</param>
  611.     /// <param name="x0">Начальное приближение (вектор)</param>
  612.     /// <param name="tol">Точность по норме градиента (по умолчанию 1e-8)</param>
  613.     /// <param name="maxIter">Максимальное число итераций (по умолчанию 20)</param>
  614.     /// <param name="eps">Шаг для численного дифференцирования (по умолчанию 1e-3)</param>
  615.     /// <param name="alpha0">Начальное значение alpha в методе Ньютона Xnext = Xcurr - alpha * H⁻¹ * ∇R(x) </param>
  616.     /// <param name="alpha_mult">Множитель на alpha, если вычисленное значение Xnext не уменьшает невязку</param>
  617.     /// <returns>Точка минимума x* (вектор)</returns>
  618.     /// 
  619.     public static double[] NewtonMinimize(Func<double[], double> R,                                              
  620.                                             double[] x0,
  621.                                             Func<double[], bool>? Validator = null,
  622.                                             double tol = 1e-8,
  623.                                             int maxIter = 40,
  624.                                             double eps = 1e-3,
  625.                                             double alpha0 = 1,
  626.                                             double alpha_mult = 0.5
  627.                                             )
  628.  
  629.     {
  630.         if (R == null) throw new ArgumentNullException(nameof(R));
  631.         if (x0 == null) throw new ArgumentNullException(nameof(x0));
  632.         if (maxIter <= 0) throw new ArgumentException("maxIter must be positive.");
  633.  
  634.         // Если функция валидации текущего приближения не задана, то по умолчанию все компоненты вектора проверяются на неотрицательность
  635.         Validator ??= a => a.All(x => x>=0); 
  636.  
  637.         double R0 = R(x0);
  638.         double Ri = R0;
  639.         double alpha = alpha0;
  640.  
  641.         double[] x = new double[x0.Length];
  642.         Array.Copy(x0, x, x0.Length);
  643.  
  644.         for (int iter = 0; iter < maxIter; iter++)
  645.         {
  646.             // Находим градиент невязки ∇R(x) (вектор частных производных невязки по компонентам вектора x)
  647.             double[] grad = Gradient(R, x, eps);
  648.             // Находим Гессиан невязки H (матрица частных вторых производных невязки по компонентам вектора x)
  649.             double[][] hess = Hessian(R, x, eps);
  650.             R0 = Ri;
  651.             Ri = R(x);
  652.  
  653.             if (Math.Abs(Ri) < eps || (Math.Abs(Ri-R0) < eps * Math.Abs(R0) && iter > 0))
  654.                 return x;
  655.  
  656.             double gradNorm = Math.Sqrt(grad.Sum(g => g * g));
  657.             if (gradNorm < tol)
  658.             {                    
  659.                 return x;
  660.             }
  661.  
  662.             // Находим -∇R(x)
  663.             double[] negGrad = AddVectorsWithScale(new double[grad.Length], grad, -1);
  664.             alpha = alpha0;
  665.             try
  666.             {
  667.                 // Решаем H * d = -∇R(x), откуда d = - H⁻¹ * ∇R(x)
  668.                 double[] d = SolveLinearSystem(hess, negGrad);
  669.  
  670.                 for (int j = 0; j < maxIter; j++)
  671.                 {
  672.                     // Xnext = Xcurr - alpha * H⁻¹ * ∇R(x)
  673.                     double[] xNew = AddVectorsWithScale(x, d, alpha);
  674.  
  675.                     if (!Validator(xNew))
  676.                     {
  677.                         alpha *= alpha_mult;
  678.                         continue;
  679.                     }
  680.  
  681.                     double RNew = R(xNew);
  682.                     // На каждой итерации невязка должна уменьшаться. Если это не так, то надо уменьшать коэффициент alpha
  683.                     if (Math.Abs(RNew) < Math.Abs(Ri))
  684.                     {
  685.                         Array.Copy(xNew, x, x.Length);
  686.                         break;
  687.                     }
  688.                     else
  689.                     {
  690.                         alpha *= alpha_mult;
  691.                     }
  692.                 }
  693.             }
  694.             catch (InvalidOperationException ex)
  695.             {
  696.                 throw new InvalidOperationException($"Newton step failed at iteration {iter}: {ex.Message}");
  697.             } 
  698.  
  699.         }
  700.  
  701.         return x;
  702.     }
  703.  
  704.  
  705.  
  706.  
  707.  
  708. }
  709.  
  710. }
  711.  
Success #stdin #stdout 0.05s 21796KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/JR3XCH
language:
C# (NET 6.0)
created:
15 hours ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!