Программа на метод сопряженных градиентов

@datafile4 · Регистрация: 26.02.2013

Студворк — интернет-сервис помощи студентам

Вот код:

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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
 
namespace ConjugateGradientVisual
{
    class Program
    {
        const double eps = 0.01;
        const double N = 2;
 
        static void Main(string[] args)
        {
            List<double> X0 = new List<double>();
            X0.Add(1);
            X0.Add(1);
            List<double> result = new List<double>((ConjugateGradient(X0)));
            Console.WriteLine("{0}  {1}", result[0], result[1]);
            Console.ReadKey();
        }
 
        static double f(List<double> X)
        {
            return 2 * X[0] * X[0] + 2 * X[1] * X[1] + 2 * X[0] * X[1] + 20 * X[0] + 10 * X[1] + 10;
        }
 
        static double f(List<double> Xk, List<double> dk, double t)
        {
            double tmp1 = Xk[0] - t * dk[0];
            double tmp2 = Xk[1] - t * dk[1];
            return 2 * tmp1 * tmp1 + 2 * tmp1 * tmp1 + 2 * tmp1 * tmp2 + 20 * tmp1 + 10 * tmp2 + 10;
        }
 
        static List<double> gradient(List<double> X)
        {
            double h = 0.0001;
            double result;
            List<double> resultV = new List<double>();
            List<double> upd = new List<double>();
            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < N; j++)
                {
                    if (j == i)
                    {
                        upd.Add(X[i] + h);
                    }
                    else {
                        upd.Add(X[i]);
                    }
                }
                result = (f(upd) - f(X)) / h;
                resultV.Add(result);
                upd.Clear();
            }
            return resultV;
        }
 
        static double GoldenMethod(double a, double b, List<double> Xk, List<double> direction)
        {
            double x = a + (3 - Math.Sqrt(5)) * (b - a) / 2;
            double y = a + (Math.Sqrt(5) - 1) * (b - a) / 2;
            double fx = f(Xk, direction, x);
            double fy = f(Xk, direction, y);
            while (b - a > eps)
            {
                if (fx <= fy)
                {
                    b = y;
                    y = x;
                    fy = fx;
                    x = a + b - y;
                    fx = f(Xk, direction, x);
                }
                else
                {
                    a = x;
                    x = y;
                    fx = fy;
                    y = a + b - x;
                    fy = f(Xk, direction, y);
                }
            }
            double point = (b + a) / 2;
            return point;
        }
        static double skalar(List<double> X, List<double> Y)
        {
            double sum = 0;
            for (int i = 0; i < N; i++)
            {
                sum += X[i] * Y[i];
            }
            return sum;
        }
 
        static double norm(List<double> X)
        {
            double sum = 0;
            for (int i = 0; i < N; i++)
                sum += X[i] * X[i];
            return Math.Sqrt(sum);
        }
 
        //Polak-Ribiere method
        static double polak(List<double> fk, List<double> fkPrev)
        {
            List<double> tmp = new List<double>();
            for (int i = 0; i < N; i++)
            {
                tmp.Add(fk[i] - fkPrev[i]);
            }
            double result = skalar(fk, tmp) / skalar(fk, fk);
            return result;
        }
 
        static List<double> direction(List<double> fk, List<double> fkPrev, double bk)
        {
            List<double> result = new List<double>();
            for (int i = 0; i < N; i++)
                result.Add(fk[i] + bk * fkPrev[i]);
            return result;
        }
 
        static List<double> ConjugateGradient(List<double> X0)
        {
            List<double> gradF = new List<double>();
            double grBk;
            List<double> fkPrev = new List<double>();
            List<double> result = new List<double>();
            List<double> xNow = new List<double>(X0);
            List<double> xPrev = new List<double>();
            double t = 0;
            List<double> grDirection = new List<double>();
            int k = 0;
 
 
            do
               {
                gradF.Clear();
                gradF.AddRange(gradient(xNow));
                grDirection.Clear();
                if (k == N || k == 0)
                {
                    grDirection.AddRange(gradF);
                    k = 0;
                }
                else {
                    fkPrev.AddRange(gradient(xPrev));
                    grBk = polak(gradF, fkPrev);
                    grDirection.AddRange(direction(gradF, fkPrev, grBk));
                }
                t = GoldenMethod(-4, 4, gradF, grDirection);
                xPrev.Clear();
                xPrev.AddRange(xNow);
                xNow.Clear();
                for (int i = 0; i < N; i++)
                {
                    xNow.Add(Math.Round((xPrev[i] - t * grDirection[i]), 4, MidpointRounding.AwayFromZero));
                }
                k++;
            } while (norm(gradF) >= eps);
            return xNow;
        }
    }
}

Делал все по алгоритму, изложенному здесь: http://math.semestr.ru/optim/fletcher.php

В конце мне выдается результат NaN и мне кажется,что это из-за значений, которые выдает функция direction. Как ни бьюсь, исправить не получается.

@afront · 10.04.2016, 14:49

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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
namespace ConGrad
{
    public delegate double MyFxDelegate(int nNumVars, ref double[] fX, ref double[] fParam);
    class Program
    {
        static void Main(string[] args)
        {
            int nNumVars = 2;
            double[] fX = new double[] { 0, 0 };
            double[] fParam = new double[] { 0, 0 };
            int nIter = 0;
            int nMaxIter = 100;
            double fEpsFx = 0.0000001;
            int i;
            double fBestF;
            string sErrorMsg = "";
            CConjugateGradient1 oOpt;
            MyFxDelegate MyFx = new MyFxDelegate(Fx3);
            oOpt = new CConjugateGradient1();
 
            for (i = 0; i < nNumVars; i++)
                {
                    Console.WriteLine("X({0}) = {1}", i + 1, fX[i]);
                }
                Console.WriteLine("Function tolerance = {0}", fEpsFx);
                Console.WriteLine("Maxumum cycles = {0}", nMaxIter);
 
            Console.WriteLine("******** FINAL RESULTS *************");
            fBestF = oOpt.CalcOptim(nNumVars, ref fX, ref fParam, fEpsFx, nMaxIter, ref nIter, ref sErrorMsg, MyFx);
 
            Console.WriteLine("Optimum at");
            for (i = 0; i < nNumVars; i++)
            {
                Console.WriteLine("X({0}) = {1}", i + 1, fX[i]);
            }
            Console.WriteLine("Function value = {0}", fBestF);
            Console.WriteLine("Number of iterations = {0}", nIter);
            Console.WriteLine();
            Console.Write("Press Enter to end the program ...");
            Console.ReadLine();
        }
        static public double Fx3(int N, ref double[] X, ref double[] fParam)
        {
          return 2 * X[0] * X[0] + 2 * X[1] * X[1] + 2 * X[0] * X[1] + 20 * X[0] + 10 * X[1] + 10;
        }
        public class CConjugateGradient1
        {
            MyFxDelegate m_MyFx;
            public double MyFxEx(int nNumVars, ref double[] fX, ref double[] fParam, ref double[] fDeltaX, double fLambda)
            {
                int i;
                double[] fXX = new double[nNumVars];
 
                for (i = 0; i < nNumVars; i++)
                {
                    fXX[i] = fX[i] + fLambda * fDeltaX[i];
                }
 
                return m_MyFx(nNumVars, ref fXX, ref fParam);
            }
 
            private void GetGradients(int nNumVars, ref double[] fX, ref double[] fParam, ref double[] fDeriv, ref double fDerivNorm)
            {
                int i;
                double fXX, H, Fp, Fm;
                fDerivNorm = 0;
                for (i = 0; i < nNumVars; i++)
                {
                    fXX = fX[i];
                    H = 0.01 * (1 + Math.Abs(fXX));
                    fX[i] = fXX + H;
                    Fp = m_MyFx(nNumVars, ref fX, ref fParam);
                    fX[i] = fXX - H;
                    Fm = m_MyFx(nNumVars, ref fX, ref fParam);
                    fX[i] = fXX;
                    fDeriv[i] = (Fp - Fm) / 2 / H;
                    fDerivNorm += Math.Pow(fDeriv[i], 2);
                }
                fDerivNorm = Math.Sqrt(fDerivNorm);
            }
 
            public bool LinSearch_DirectSearch(int nNumVars, ref double[] fX, ref double[] fParam, ref double fLambda, ref double[] fDeltaX, double InitStep, double MinStep)
            {
                double F1, F2;
 
                F1 = MyFxEx(nNumVars, ref fX, ref fParam, ref fDeltaX, fLambda);
 
                do
                {
                    F2 = MyFxEx(nNumVars, ref fX, ref fParam, ref fDeltaX, fLambda + InitStep);
                    if (F2 < F1)
                    {
                        F1 = F2;
                        fLambda += InitStep;
                    }
                    else
                    {
                        F2 = MyFxEx(nNumVars, ref fX, ref fParam, ref fDeltaX, fLambda - InitStep);
                        if (F2 < F1)
                        {
                            F1 = F2;
                            fLambda -= InitStep;
                        }
                        else
                        {
                            // reduce search step size
                            InitStep /= 10;
                        }
                    }
                } while (!(InitStep < MinStep));
 
                return true;
 
            }
 
 
            public double CalcOptim(int nNumVars, ref double[] fX, ref double[] fParam, double fEpsFx, int nMaxIter, ref int nIter, ref string sErrorMsg, MyFxDelegate MyFx)
            {
 
                int i;
                double[] fDeriv = new double[nNumVars];
                double[] fDerivOld = new double[nNumVars];
                double F, fDFNormOld, fLambda, fLastF, fDFNorm = 0;
 
                m_MyFx = MyFx;
 
                // calculate and function value at initial point
                fLastF = MyFx(nNumVars, ref fX, ref fParam);
 
                GetGradients(nNumVars, ref fX, ref fParam, ref fDeriv, ref fDFNorm);
 
                fLambda = 0.1;
                if (LinSearch_DirectSearch(nNumVars, ref fX, ref fParam, ref fLambda, ref fDeriv, 0.1, 0.000001))
                {
                    for (i = 0; i < nNumVars; i++)
                    {
                        fX[i] += fLambda * fDeriv[i];
                    }
                }
                else
                {
                    sErrorMsg = "Failed linear search";
                    return fLastF;
                }
 
                nIter = 1;
                do
                {
                    nIter++;
                    if (nIter > nMaxIter)
                    {
                        sErrorMsg = "Reached maximum iterations limit";
                        break;
                    }
                    fDFNormOld = fDFNorm;
                    for (i = 0; i < nNumVars; i++)
                    {
                        fDerivOld[i] = fDeriv[i]; // save old gradient
                    }
                    GetGradients(nNumVars, ref fX, ref fParam, ref fDeriv, ref fDFNorm);
                    for (i = 0; i < nNumVars; i++)
                    {
                        fDeriv[i] = Math.Pow((fDFNorm / fDFNormOld), 2) * fDerivOld[i] - fDeriv[i];
                    }
                    if (fDFNorm <= fEpsFx)
                    {
                        sErrorMsg = "Gradient norm meets convergence criteria";
                        break;
                    }
                    fLambda = 0;
                    if (LinSearch_DirectSearch(nNumVars, ref fX, ref fParam, ref fLambda, ref fDeriv, 0.1, 0.000001))
                    {
                        for (i = 0; i < nNumVars; i++)
                        {
                            fX[i] += fLambda * fDeriv[i];
                        }
                        F = MyFx(nNumVars, ref fX, ref fParam);
                        if (Math.Abs(F - fLastF) < fEpsFx)
                        {
                            sErrorMsg = "Successive function values meet convergence criteria";
                            break;
                        }
                        else
                        {
                            fLastF = F;
                        }
                    }
                    else
                    {
                        sErrorMsg = "Failed linear search";
                        break;
                    }
                } while (true);
 
                return fLastF;
            }
        }
    }
}

Новые блоги и статьи Все статьи Все блоги /
Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Programma_Boinc 28.12.2025 Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Налог на собак: https:/ / ********/ gallery/ V06K53e Финансовый отчет в Excel: https:/ / ********/ gallery/ bKBkQFf Пост отсюда. . .	Кто-нибудь знает, где можно бесплатно получить настольный компьютер или ноутбук? США. Programma_Boinc 26.12.2025 Нашел на реддите интересную статью под названием Anyone know where to get a free Desktop or Laptop? Ниже её машинный перевод. После долгих разбирательств я наконец-то вернула себе. . .	Thinkpad X220 Tablet — это лучший бюджетный ноутбук для учёбы, точка. Programma_Boinc 23.12.2025 Рецензия / Мнение/ Перевод Нашел на реддите интересную статью под названием The Thinkpad X220 Tablet is the best budget school laptop period . Ниже её машинный перевод. Thinkpad X220 Tablet —. . .	PhpStorm 2025.3: WSL Terminal всегда стартует в ~ and_y87 14.12.2025 PhpStorm 2025. 3: WSL Terminal всегда стартует в ~ (home), игнорируя директорию проекта Симптом: После обновления до PhpStorm 2025. 3 встроенный терминал WSL открывается в домашней директории. . .	Как объединить две одинаковые БД Access с разными данными VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.
Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .