Задача
Найти минимальное значение функции F(x)=f(x)+exp(-x*x), где f(x) задана дифференциальным уравнением
Первый взгяд
Нам надо найти минимальное значение функции, в задании ничего не сказано про отрезок, на котором нужно искать. На ум приходит найти точку минимума, если она есть. Это и будет минимальное значение функции.
Поиск экстремума
Рассмотрим функцию аналитически: найдем производную и приравняем ее нулю. Тем самым, найдем точки экстремума. Так как одна часть функции уже задана дифф уравнением, найдем производную от второй части — экспоненты и приравняем все нулю, затем найдем корни уравнения:
Чтобы определить, возрастает или убывает функция на интервале, необходимо в производную подставить любое число из этого интервала. Получаем следующее
Если знак производной меняет свой знак с «+» на «-«, значит, точка максимума. Дела… Получается, что если Х<0.39 — функция возрастает и минимальное значение будет при Х=0, и если Х>0.39 — то ближе к правой границе отрезка, на котором мы рассматриваем функцию.
Алгоритм
Для поиска минимального значения будем использовать самый простой метод — дробление шага. Сложность данной задачи состоит в том, что конечной функции, для которой можно искать мин. значение, у нас нет. У нас есть часть конечной функции, заданная как дифф. уравнение, поэтому мы должны для каждого X решить дифф. уравнение, и для этого X найти значение конечной функции, потом сравнить с предыдущим рассчитанным значением и выполнить либо дробление шага, либо шаг.
Для решения дифф. уравнения будем использовать метод Рунге-Кутта 4 порядка. Почитать про метод тут. Так как решение дифф. уравнения зависит от начального условия, поэтому каждый раз решение дифф. уравнения для каждого X будем считать заново, с начального условия f(0)=0.
Программа
Программа решения дифф. уравнения методом Рунге-Кутты 4 порядка
| 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 | // пересчет правых частей уравнений // n - размерность задачи // х - текущее значение x // y - текущее значение левых частей дифф. уравнений y // f - значение правой части функции f(x,y) void prav(int n, double x, double *y, double *f) { f[1] = 2.0/sqrt(2.0*pi)*exp(-x*x); // правая часть дифф. уравнения } // метод Рунге-Кутта 4 степени // n - размерность задачи // x - текущий x до шага // y - значение левой части уравнения void RKM(int n, double &x, double *y, double h) { double *f = new double[n + 1]; double *y1 = new double[n + 1]; double *fk = new double[5, n + 1]; prav(n, x, y, f); //Расчет правых частей for (int i = 1; i <= n; i++) // в точке X { fk[1, i] = h * f[i]; // значение K1 y1[i] = y[i] + 0.5 * fk[1, i]; } x = x + 0.5 * h; prav(n, x, y1, f); for (int i = 1; i <= n; i++) // в точке X + 0.5 * h { fk[2, i] = h * f[i]; // значение k2 y1[i] = y[i] + 0.5 * fk[2, i]; } prav(n, x, y1, f); for (int i = 1; i <= n; i++) // в точке X + 0.5 * h { fk[3, i] = h * f[i]; // значение k3 y1[i] = y[i] + fk[3, i]; } x = x + 0.5 * h; prav(n, x, y1, f); for (int i = 1; i <= n; i++) // в точке X + h { fk[4, i] = h * f[i]; // значение k4 // решение в точке x + h y[i] = y[i] + 0.1666667 * (fk[1, i] + 2 * fk[2, i] + 2 * fk[3, i] + fk[4, i]); } return; } |
Для правильного решения дифф. уравнения нужно учитывать начальное условие, поэтому каждый новое значение x считаем заново с начального условия
| 1 2 3 4 5 6 7 8 9 10 | double rkm_x(double x, double h) // метод рунге-кутты для значения x с учетом начальное условия { double xx=0; // начальное условие x=0 double *y = new double[2]; y[1]=0; // начальное условия y(0)=0 while(xx<=x-h) { // от хх=0 до текущего х RKM(1,xx,y,h); // считаем методом рунге-кутты в точке x } return y[1]; } |
После всех этих функций мы имеем решение дифф. уравнения в точку Х, найдем значение конечной функции F с помощью функции
| 1 2 3 4 5 6 7 | // f - конечная функция // f1 - посчитанный дифф. уравнение в точке х // х - значение в точку Х void func(double &f, double f1, double x) { f = f1+exp(-x*x); } |
Основная программа, которая реализует метод дробления шага
| 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 | int _tmain(int argc, _TCHAR* argv[]) { double a,b; // границы отрезка, на котором ищем минимальное значение double h; // шаг cout << "Vvedite polozhitelnie granici otrezka [A,B], na kotorom nyzhno najti min znachenie funkcii:" << endl;// вывод на экран cout << "Levaya granica A: "; // вывод на экран cin >> a; // ввод с клавиатуры cout << "Pravaya granica B: "; // вывод на экран cin >> b; // ввод с клавиатуры cout << "vvedite shag h: "; // вывод на экран cin >> h; // ввод с клавиатуры cout << endl; // вывод на экран с новой строку int n=1; // размерность задачи (у нас одна функция и одно уравнение дифф.) double *y = new double[2]; // в массиве y[0] не учитывается double f=0, fpred=0; // значение функции в точке и предыдущее значение double eps = 0.001; // точность поиска минимального значение (можно задавать с клавы) //fstream fileA("y(x).txt", ios_base::in|ios_base::out|ios_base::trunc); //fileA << "x y" << endl; double x=a; // начинаем с самой левой границы double xpred=0; // предыд. шаг func(f,0,0); // считаем в точке 0 cout << "f(0) = " << f << endl; // вывод на экран cout << "Na4alnij x: " << x << endl; cout << "Na4alnij h: " << h << endl; double f1 = rkm_x(x,h); // в точке x+h func(f,f1,x); // считаем в точку x+h fpred=f; xpred = x; cout << "f(" << x << ") = " << f << endl; // вывод на экран cout << "===========================" << endl; while(true) { while (x+h>=b) // проверка выхода за правую границу с учетом сделанного шага { h /= 2; //делим шаг пополам и опять пробуем выход за границу } // x += h; // делаем шаг double f1 = rkm_x(x,h); // в точке x func(f,f1,x); // считаем в точку x cout << "tek. x: " << x << endl; cout << "tek. h: " << h << endl; cout << "pred. f(" << xpred << "): " << fpred << endl; cout << "tek. f(" << x << "): " << f << endl; if (fpred<f) cout << "f(" << xpred << ") < f(" << x << ") - delim h popolam" << endl; else cout << "f(" << xpred << ") > f(" << x << ")" << endl; cout << "Ysl. abs(f-fpred): " << abs(f-fpred) << endl; if (abs(f-fpred)<eps || abs(x-xpred)<eps) //пока значение текущего и предыдущего не выходит за точность break; if (fpred<f) { //если текущее значение функции больше предыдущего x -= h; // то вычитаем шаг if (x <= a) x=a; // проверка на выход за границы h /= 2.0; // делим на шаг } else { fpred = f; //сохраняем предыдущее значение функции f xpred = x; } cout << "Tek. x v konce cikla: " << x << endl; // вывод на экран cout << "Tek. h v konce cikla: " << h << endl; // вывод на экран if (fpred<f) cout << "f(" << xpred << ") v konce cikla: " << fpred << endl; // вывод на экран else cout << "f(" << x << ") v konce cikla: " << f << endl; // вывод на экран cout << endl; // новая строчка system("PAUSE"); // паузка жмем ентер } cout << "===========================" << endl; cout << "Min F(" << xpred << ") = " << fpred << endl; // вывод значения функции в точке x system("PAUSE"); // жмем ентер return 0; // если да, то выход } |
Скриншоты
При X<0.39
При X>0.39
