Форум программистов, компьютерный форум, киберфорум
Наши страницы
С++ для начинающих
Войти
Регистрация
Восстановить пароль
 
 
Рейтинг 4.86/7: Рейтинг темы: голосов - 7, средняя оценка - 4.86
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
1

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

16.09.2017, 21:29. Просмотров 1303. Ответов 22

Доброго дня суток! Очень нужна ваша помощь! Необходимо минимизировать методом Ньютона функции Розенброка и Химмельблау , т.е.
Минимизация двумерным методом Ньютона Функции Химмельблау и Розенброка


Написал код,но для этих двух функций он не работает. Как его оптимизировать под эти функции? спасибо заранее4

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#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;
0
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
16.09.2017, 21:29
Ответы с готовыми решениями:

Минимизация Методом Ньютона уравнения с двумя переменными
Всем добрый день... я тут новичок...очень нужна ваша помощь. Программу для уравнения где одна...

Нахождения максимума функции методом Ньютона
Мне надо на с++ написать код для нахождения максимума функции методом ньютона. Код бы я написать...

Как найти максимум функции методом Ньютона?
кто может помочь? Я написал программу для поиска корней а вот для поиска максимума на отрезке не...

Ошибка в программе интерполяция функции методом Ньютона
выдает ошибку, не пойму в чем проблема Unit1.cpp(12): E2313 Constant expression required ...

Интерполяция функции одной переменной методом Ньютона (Pascal -> C++)
program interpol; uses crt,graph; const MAXCOUNT=30; type

22
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
17.09.2017, 11:34 2
Ты что-то непонятное написал какие-то циклы страшные, там один цикл - итерации методом Ньютона, матрица 2 на 2 - ее легко инвертировать, функции хорошие, метод сходится с любого приближения, вот моя реализация. Только вывожу результат в файлы, а не на консоль. Может якобиан градиента проверить стоит ибо в спешке считал, но вроде все нормально сходится.
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#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();
}
Добавлено через 3 часа 37 минут
Поправил опечатку.
0
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
17.09.2017, 15:33  [ТС] 3
не компилится к сожалению...ошибок море выдал (
0
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
17.09.2017, 15:43 4
IcebergLife, Что за ошибки? Тут особо никаких ошибок быть не может.
1
17.09.2017, 15:43
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
17.09.2017, 15:47  [ТС] 5
ну вот пример :
23 cтрока ostream не является членом std, os необъявленный идентификатор, и тд...
и так в многих строчка в основном что например os необъявленный идентификатор и endl,ofstream и ostream не является членом std
0
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
17.09.2017, 15:53 6
Если VS, добавь #include "stdafx.h" первой строчкой, может поможет
1
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
17.09.2017, 16:08  [ТС] 7
да я два инклуда добавил не помогает
#include "stdafx.h"
#include <stdio.h>

Добавлено через 1 минуту
и добавил using namespace std;

Добавлено через 1 минуту
а все, спасибо. добавил в самое начало - заработало.

Добавлено через 2 минуты
еще вопрос. если мне нужно привязать к критерию остановы e<=0.1,0.01,0,001 , как это в программу внедрить?

Добавлено через 9 минут
еще вопрос. если мне нужно привязать к критерию остановы e<=0.1,0.01,0,001 , как это в программу внедрить?
0
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
17.09.2017, 16:14 8
Всмысле тебе нужно задавать значение e?
Тогда функция будет такой:
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// Поиск экстремального значения методом Ньютона
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)
0
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
17.09.2017, 16:25  [ТС] 9
спасибо еще раз! большущее!
0
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
21.09.2017, 11:01  [ТС] 10
Дорогой друг,привет еще раз. Опять нужна твоя помощь. очень. по этому коду.

можно ли кроме точек x1 и x2,которые есть рассмотреть другие варианты? почему именно эти выбираются?
когда я считал функцию розенброка параметр точности е не влияет вообще ни на что почему то...
для функции Химмельблау ситуация аналогична и у этой функции 4 минимума и надо как то остальные рассмотреть...

спасибо заранее
0
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
21.09.2017, 11:08 11
Цитата Сообщение от IcebergLife Посмотреть сообщение
почему именно эти выбираются?
Скорее всего, они выбираются исходя из начального приближения. Вот здесь:
C++
1
Vector2 rx(10.0, 10.0); // Начальное значение
Находится ближайшее к заданному.
Теоретически тебе нужно приближенно найти точку вблизи нужного минимума и передать методу, он уточнит ее значение до заданного e. Как найти начальное значение - я не знаю.
0
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
21.09.2017, 11:12  [ТС] 12
мб e не влияет потому что вообще элементы матрицы неправильно посчитаны? или исключено?

а как помимо функции розенброка посчитать это все для простой квадратичной функции (x1)^2 + (x2)^2 ?
0
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
21.09.2017, 11:33 13
Цитата Сообщение от IcebergLife Посмотреть сообщение
e не влияет
А на что он должен влиять? Это же просто точность с которой ты значения считаешь или что-то другое?

Добавлено через 10 минут
IcebergLife, вот смотри: я запускаю с другого начального приближения для химмельблау:
C++
1
2
Vector2 hx(10.0, -10.0);
newton_minimization(hfunc, hx,  iteration_count, 1e-8, os);
Получаю другую точку:
Iteration: 9
(x1,x2) = 3.58443 -1.84813
Functions value: 1.57772e-29

Т.е. метод Ньютона он только уточняет значения, начальное приближение нужно находить другим способом, можно просто перебирать из некоторой области пробные точки.

Добавлено через 7 минут
Я понял! Тебе нужно брать начальные приблоижения из разных координатных четвертей (x,x) (-x,x), (x,-x), (-x,-x) и ты получишь все точки! И самое главное исправь в функции newton_minimization строчку if(f.gradient().x1 <= e && f.gradient().x2 <= e) { на строчку if(fabs(f.gradient().x1) <= e && fabs(f.gradient().x2) <= e) {, это моя ошибка
0
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
21.09.2017, 11:44  [ТС] 14
woldemas, да, e это точность и как я понимаю она должна влиять.... сейчас попробую пересчитать все.

и вопрос. можешь еще добавить в код для функции (x1)^2 + (x2)^2 такой же расчет? сделал бы сам да до сих пор не въезжаю как считать якобиана матрицу
0
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
21.09.2017, 12:42 15
Цитата Сообщение от IcebergLife Посмотреть сообщение
можешь еще добавить в код для функции (x1)^2 + (x2)^2 такой же расчет?
А что там минимизировать, там (0,0) - минимум получается.
Ну так и быть. Добавляешь этот класс:
C++
1
2
3
4
5
6
7
8
9
10
11
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;
}
И вызов метода такой же, например:
C++
1
2
3
    SquaredSum func;
    Vector2 x(10.0, 5.0); 
    newton_minimization(func, x, iteration_count, 1e-8, std::cout);
0
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
22.09.2017, 11:39  [ТС] 16
Мне сказали что надо выбирать точки так...

Для Розенброка: одна вблизи минимума, вторая на противоположном крае дна оврага, по одной на склонах оврагов, одна удаленная от оврага точка,

Для Химмельблау: по одной вблизи с каждым из минимумов, между каждой парой минимумов, в центре четырех минимумов, удаленная точка,

2. Не понимаю, что определяет точность е, на что она влияет, на координаты точки минимума, на значение ф-ции должна быть хочется заданная точность...т.е. Е=0.1 точно изменяется с десятков 0.001 с сотен и тд...
0
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
22.09.2017, 12:20  [ТС] 17
Вот во вложении скрин...ешка (точность ни на что не повлияла)..
0
Миниатюры
Минимизация двумерным методом Ньютона Функции Химмельблау и Розенброка  
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
22.09.2017, 12:58 18
Цитата Сообщение от IcebergLife Посмотреть сообщение
точность ни на что не повлияла
Она влияет на найденное значение: при e=0.0001 x1 = 1, x2 = 1; при e=0.1 x1=0.998, x2=0.994. Т.е. значения x определяются менее точно и все, больше ни на что она и не должна влиять.
0
IcebergLife
0 / 0 / 0
Регистрация: 01.06.2017
Сообщений: 18
22.09.2017, 13:04  [ТС] 19
Оно посути должно влиять например там про 0.1 3.100000 а при 0.001 3.1230000 ( это я пример привёл) мб это е как то по другому в коде нужно сделать . Примерно понимаю что но не понимаю как в коде реализовать
0
woldemas
451 / 321 / 176
Регистрация: 06.09.2013
Сообщений: 974
22.09.2017, 13:12 20
Цитата Сообщение от IcebergLife Посмотреть сообщение
мб это е как то по другому в коде нужно сделать
Может быть, попробуй вот так:
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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;
}
1
22.09.2017, 13:12
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
22.09.2017, 13:12

Минимизация функции методом ломаных
Доброго времени суток! Помогите, пожалуйста, с программой. Мне нужно реализовать метод ломаных....

Минимизация функции методом ломаных
Здравствуйте! Помогите, пожалуйста, реализовать метод ломаных нахождения минимума...

Минимизация функции методом золотого сечения
Здравствуйте, помогите с лабораторной. Нужно написать код, исходя из следующей блок схемы:...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
20
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Рейтинг@Mail.ru