Форум программистов [Powered by Invision Power Board]

Цитата(MrNull @ 23.12.2007, 15:42)

Делюсь с остальными:

Программа поиска безусловного минимума функции f(x1, x2) = 2 * (x1 - 3)^2 + (x2 - 2)^2

    Алгоритм расчета: методом Ньютона решается система
            |F1(x1, x2) = df/dx1 = 0|
    (1)        |                        |
            |F2(x1, x2) = df/dx2 = 0|,
    состоящая из равенства нулю частных производных исходной функции
    (необходимое условие экстремума).

    Алгоритм метода Ньютона:
        Шаг 1. Задается начальное приближение для x1, x2,...xn (в нашем случае n = 2)
        Шаг 2. Вычисляется невязка fvec (вектор правых частей системы (1))
               Проверяется сходимость (сумма модулей координат fvec меньше некоторого малого eps > 0)
               Если достигнута сходимость, цикл завершается. Иначе переход к следующему шагу.
        Шаг 3. Вычисляются элементы матрицы Якобиана системы (1):

                     |dF1(x1, x2)/dx1    dF1(x1, x2)/dx2|
               J  =  |                                   |
                     |dF2(x1, x2)/dx1    dF2(x1, x2)/dx2|

        Шаг 4. Методом Гаусса решается система J * (dx) = - fvec, где dx = (dx1, dx2) - вектор смещений
        Шаг 5. Рассчитывается следующее приближение переменных:
               x1 = x1 + dx1,    x2 = x2 + dx2
        Шаг 6. Переход к следующей итерации (Шаг 1.) если не превышено максимальное число итераций

    Метод Гаусса решения линейных систем:
        Элементарными преобразованиями над строками
        матрицы она приводится к верхнетреугольному виду (прямой ход метода Гаусса);
        полученная система с треугольной матрицей решается в явном виде
        методом последовательного исключения неизвестных в порядке dxn,...,dx1
        (обратный ход метода Гаусса)

Код


int main()
{        
    const int n = 2;            //размерность вектора    
    double x[n];                //вектор координат исходного минимума
    const int ntrial = 10;        //максимально число шагов для метода Ньютона
    const double eps = 1.e-6;    //точность (см. Шаг 2)
    int iter;                    //номер итерационного шага в методе Ньютона
    int i, j, k;                //переменные циклов
    double tmp;                    //временная переменная
    double fvec[n];                //вектор для хранения частных производных (см. Шаг 2)
    double p[n];                //вектор для хранения вектора смещений (см. Шаг 4)
    double fjac[n][n];            //матрица Якобиана (см. Шаг 3)
    
    //Создаем файл для вывода результатов расчета (указать полный путь)
    const char *logname = "C:\\Log.txt";
    ofstream os(logname);
    if (!os) {
        cout << "Cannot create log file " << logname << "\n";
        return 0;
    }

//Для вывода на экран
#define stream cout
//Для вывода в файл
//#define stream os

    //Шаг 1
    //Можно задать любое начальное приближение
    x[0] = 123;
    x[1] = 321;

    //Итерационный процесс метода Ньютона
    for (iter = 1; iter <= ntrial; iter++) {
        stream << "\n========================== " << iter << " iteration";
        
        stream << "\n\nInfo: Current vector x:";
        for (i = 0; i < n; i++)
            stream << "\nx[" << i << "] = " << x[i];

        //Шаг 2  

        fvec[0] = 4 * (x[0] - 3);    //df/dx1 = 4 * (x1 - 3)
        fvec[1] = 2 * (x[1] - 2);    //df/dx2 = 2 * (x2 - 2)
                 
        tmp = 0.; 
        for (i = 0; i < n; i++) { 
            tmp += fabs(fvec[i]);
        }
        if (tmp <= eps) 
            break;
        
        //Шаг 3

        fjac[0][0] = 4;        //dF1(x1, x2)/dx1 = 4        
        fjac[0][1] = 0.;    //dF1(x1, x2)/dx1 = 0
        fjac[1][0] = 0.;    //dF1(x1, x2)/dx1 = 0    
        fjac[1][1] = 2.;    //dF1(x1, x2)/dx2 = 2

        stream << "\n\nInfo: Current Jacobean matrix at vector x:\n";
        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                stream << fjac[i][j] << "\t";
            }            
            stream << "\n";
        }

        //Шаг 4

        //Правая часть линейной системы 
        for (i = 0; i < n; i++) 
            p[i] = - fvec[i]; 

        //Прямой ход метода Гаусса
        for (i = 0; i < n - 1; i++) {
            tmp = fjac[i][i];
            
            if (tmp == 0.) {
                stream << "Error: Zero element at main diagonal in Jacobean matrix\n";
                return 0;
            }
            
            // рабочий элемент корректен - обрабатываем все нижлежащие строки
            for (j = i + 1; j < n; j++) {
                p[j] -= p[i] / tmp;
                for (k = n - 1; k >= i; k--) {
                    fjac[j][k] -= fjac[i][k] * fjac[j][i] / tmp;
                }
            }            
        }

        //Нахождение вектора смещений (обратный ход метода Гаусса)
        for (i = n - 1; i >= 0; i--) {
            tmp = 0;
            for (j = i + 1; j < n; j++) {
                tmp += fjac[i][j] * p[j];
            }
            p[i] = (p[i] - tmp) / fjac[i][i];
        }
        
        //Шаг 5
        for (i = 0; i < n; i++) { 
            x[i] += p[i];
        }
        
        //Шаг 6
    }

    if (iter < ntrial) {
        stream << "\n\n========================== The minimum location\n";
        stream << "(" << x[0] << ", " << x[1] << ")\n";
        return 0;
    }    
    stream << "\nError: Maximum number of iterations exceed\n";
    return 0;
}