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

#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();
}