Минимизация двумерным методом Ньютона Функции Химмельблау и Розенброка

#include "stdafx.h"
#include <math.h>
#include <conio.h>
#include <stdio.h>
#include <iostream>
#include <cmath>
#include <string.h>
 
using namespace std;
 
 
int main()
 
{
    const int n = 2;            
    double x[n];                
    const int shag = 40;        
    int iter;                    
 
       int i, j, k;                
    double tmp;                    
    double fvec[n];                
    double p[n];                
    double fmatrix[n][n];            
 
    double e;               std::cout << "Vvedite e" << std::endl;
    std::cin >> e;
#define stream cout
 
    x[0] = 10;
    x[1] = 30;
 
 
for (iter = 1; iter <= shag; iter++) {
        stream << "\n " << iter << " iteraciya";
 
stream << "\n\nInfo: Vektor x:";
for (i = 0; i < n; i++)
stream << "\nx[" << i << "] = " << x[i]; 
 
 
 
 
fvec[0] = 2 * x[0];    
fvec[1] = 2 * x[1];    
 
tmp = 0.;
for (i = 0; i < n; i++) {
tmp += fabs(fvec[i]);
        }
if (tmp <= e)
break;
 
 
fmatrix[0][0] = 2;        
fmatrix[0][1] = 0;    
fmatrix[1][0] = 0;    
fmatrix[1][1] = 2;    
 
stream << "\n\nInfo: Vectornaya matrica x:\n";
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
stream << fmatrix[i][j] << "\t";
            }
stream << "\n";
        }
 
        for (i = 0; i < n; i++)
            p[i] = -fvec[i];
            for (i = 0; i < n - 1; i++) {
            tmp = fmatrix[i][i];
 
            if (tmp == 0.) {
                stream << "Error: Nulevoy element v matritse\n";
    return 0;
            }
 
 
for (j = i + 1; j < n; j++) {
p[j] -= p[i] / tmp;
for (k = n - 1; k >= i; k--) {
                    fmatrix[j][k] -= fmatrix[i][k] * fmatrix[j][i] / tmp;
                }
            }
        }   
for (i = n - 1; i >= 0; i--) {
tmp = 0;
for (j = i + 1; j < n; j++) {
tmp += fmatrix[i][j] * p[j];
            }
p[i] = (p[i] - tmp) / fmatrix[i][i];
        }
 
 
for (i = 0; i < n; i++) {
x[i] += p[i];
        }
    }
if (iter < shag) {
        stream << "\n\n Minimum \n";
        stream << "(" << x[0] << ", " << x[1] << ")\n";
        return 0;
    }
    stream << "\nError: Maksimalnoe chislo iteracii\n";
    return 0;

#include <iostream>
#include <fstream>
#include <cmath>
 
struct Vector2 {
    double x1;
    double x2;
    Vector2() : x1(0.0), x2(0.0) {}
    Vector2(double _x1, double _x2): x1(_x1), x2(_x2) {}
    Vector2 & operator += (const Vector2 &other) {
        x1 += other.x1;
        x2 += other.x2;
        return *this;
    }
    Vector2 operator-() const {   return Vector2(-x1, -x2);  }    
};
 
std::ostream & operator<<(std::ostream & os, const Vector2 &x) { return os << x.x1 << " " << x.x2; }
 
class Matrix2x2 {
    double _items[2][2];
public:
    Matrix2x2() {}
    Matrix2x2(double m11, double m12, double m21, double m22) {
        _items[0][0] = m11; _items[0][1] = m12;
        _items[1][0] = m21; _items[1][1] = m22;
    }
    Vector2 operator*(const Vector2 &v) const {
        return Vector2(_items[0][0] * v.x1 + _items[0][1] * v.x2,
                       _items[1][0] * v.x1 + _items[1][1] * v.x2);
    }
    Matrix2x2 inverse() const {
        double det = _items[0][0] * _items[1][1] - _items[0][1] * _items[1][0];
        if(det == 0.0) throw std::exception();
        return Matrix2x2(_items[1][1] / det, -_items[0][1] / det,
                       -_items[1][0] / det, _items[0][0] / det);
    }
    double* operator[](size_t i) {return _items[i]; }
    const double* operator[](size_t i) const {return _items[i]; }
};
 
// Базовый класс для функции от 2-х переменных
class Func2 {
protected:
    double _value;
    Vector2 _gradient;
    Matrix2x2 _jacobian;
public:
    // Вычислить значение в данной точке
    virtual void evaluate(const Vector2 &x) = 0;
    // Значение функции
    double value() const {return _value;}
    // Градиент
    const Vector2 & gradient() const {return _gradient; }
    // Якобиан градиента
    const Matrix2x2 & jacobian() const { return _jacobian; }
 
};
 
// Функция Розенброка
class RozenbrokFunc : public Func2 {
    virtual void evaluate(const Vector2 &x);
};
 
void RozenbrokFunc::evaluate(const Vector2 &x) {
        double a = x.x1 - x.x2 * x.x2;
        double b = 1.0 - x.x2;
        _value = 100.0 * a * a + b * b;
        _gradient.x1 = 200.0 * a;
        _gradient.x2 = -800.0 * a * x.x2 - 2.0 * b;
 
        _jacobian[0][0] = 200.0;
        _jacobian[0][1] = -400.0 * x.x2;
        _jacobian[1][0] = -800.0 * x.x2;
        _jacobian[1][1] = -800.0 * (x.x1 - 3.0 * x.x2 * x.x2) + 2;
}
 
// Функция Химмельблау
class HimmelblauFunc : public Func2 {
    virtual void evaluate(const Vector2 &x);
};
 
void HimmelblauFunc::evaluate(const Vector2 &x) {
    double a = x.x1 * x.x1 + x.x2 - 11;
    double b = x.x1 + x.x2 * x.x2 - 7;
    _value = a * a + b * b;
 
    _gradient.x1 = 4.0 * a * x.x1 + 2.0 * b;
    _gradient.x2 = 2.0 * a + 4.0 * b * x.x2;
 
    _jacobian[0][0] = 2.0 + 4.0 * (2.0 * x.x1 * x.x1 + a);
    _jacobian[1][0] = _jacobian[0][1] = 4.0 * (x.x1 + x.x2);
    _jacobian[1][1] = 2.0 + 4.0 * (2.0 * x.x2 * x.x2 + b);
}
 
// Поиск экстремального значения методом Ньютона
bool newton_minimization(Func2 &f, Vector2 &x, size_t iteration_count, std::ostream & os) {
    double threshold = 0.1;
    for(size_t i = 0; i < iteration_count; i++) {
        f.evaluate(x);
        Matrix2x2 inv_j = f.jacobian().inverse();
        Vector2 xstep = - (inv_j * f.gradient());
        os << "Iteration: " << i << std::endl;
        os << "(x1,x2) = " << x << std::endl;
        os << "Functions value: " << f.value() << std::endl;
        os << "Gradient: " << f.gradient() << std::endl;
        os << "Step: " << xstep << std::endl;
        if(f.gradient().x1 <= threshold && f.gradient().x2 <= threshold) {
            os << "Solution founded" << std::endl;
            return true;
        }
        x += xstep;
        threshold *= 0.1;
        os << std::endl;
    }
    return false;
}
 
 
int main()
{
    Vector2 rx(10.0, 10.0); // Начальное значение
    size_t iteration_count = 15;
    RozenbrokFunc rfunc;
    std::ofstream os;
    os.open("rozenbrok.txt");
    newton_minimization(rfunc, rx,  iteration_count, os);
    os.close();
    os.open("himelblau.txt");
    HimmelblauFunc hfunc;
    Vector2 hx(10.0, 10.0);
    newton_minimization(hfunc, hx,  iteration_count, os);
    os.close();
}

// Поиск экстремального значения методом Ньютона
bool newton_minimization(Func2 &f, Vector2 &x, size_t iteration_count, double e, std::ostream & os) {    
    for(size_t i = 0; i < iteration_count; i++) {
        f.evaluate(x);
        Matrix2x2 inv_j = f.jacobian().inverse();
        Vector2 xstep = - (inv_j * f.gradient());
        os << "Iteration: " << i << std::endl;
        os << "(x1,x2) = " << x << std::endl;
        os << "Functions value: " << f.value() << std::endl;
        os << "Gradient: " << f.gradient() << std::endl;
        os << "Step: " << xstep << std::endl;
        if(f.gradient().x1 <= e && f.gradient().x2 <= e) {
            os << "Solution founded" << std::endl;
            return true;
        }
        x += xstep;        
        os << std::endl;
    }
    return false;
}
 
// А в вызов добавь значение, например 0.1
newton_minimization(rfunc, rx,  iteration_count, 0.1, os)

Vector2 rx(10.0, 10.0); // Начальное значение

Vector2 hx(10.0, -10.0);
newton_minimization(hfunc, hx,  iteration_count, 1e-8, os);

class SquaredSum : public Func2 {
    virtual void evaluate(const Vector2 &x);
};
 
void SquaredSum::evaluate(const Vector2 &x) {
    _value = x.x1 * x.x1 + x.x2 * x.x2;
    _gradient.x1 = 2.0 * x.x1;
    _gradient.x2 = 2.0 * x.x2;
    _jacobian[0][0] = 2.0; _jacobian[0][1] = 0.0;
    _jacobian[1][0] = 0.0; _jacobian[1][1] = 2.0;
}

    SquaredSum func;
    Vector2 x(10.0, 5.0); 
    newton_minimization(func, x, iteration_count, 1e-8, std::cout);

bool newton_minimization(Func2 &f, Vector2 &x, size_t iteration_count, double e, std::ostream & os) {
    for(size_t i = 0; i < iteration_count; i++) {
        f.evaluate(x);
        Matrix2x2 inv_j = f.jacobian().inverse();
        Vector2 xstep = - (inv_j * f.gradient());
        os << "Iteration: " << i << std::endl;
        os << "(x1,x2) = " << x << std::endl;
        os << "Functions value: " << f.value() << std::endl;
        os << "Gradient: " << f.gradient() << std::endl;
        os << "Step: " << xstep << std::endl;
        if(fabs(xstep.x1) <= e && fabs(xstep.x2) <= e) {
            os << "Solution founded" << std::endl;
            return true;
        }
        x += xstep;
        os << std::endl;
    }
    return false;
}