5. Многомерная оптимизация
5.1 Метод Ньютона
5.2 Нахождение экстремума для заданной функции
Экстремум был найден с помощью программы, написанной на языке Си. Текст программы приведён в приложении А.
Результаты вычислений:
Рисунок 5.1 — Нахождение экстремума с помощью метода Ньютона
6. Методы обработки экспериментальных данных
6.1 Общие сведения
Пусть в процессе эксперимента путем измерений получена таблица некоторой зависимости f :
Рисунок 6.1 – Исходные данные к задаче обработки экспериментальных данных в общем виде
Требуется найти формулу, выражающую эту зависимость аналитически. Сформулируем задачу таким образом, чтобы учитывался характер исходной функции: найти функцию F заданного вида y=F(x), принимающую в точках , ,…, , значения близкие к табличным значениям , ,…, .
На практике вид приближающей функции F определяется следующим образом. По заданной таблице значений f (x) строится ее точечный график, а потом проводится гладкая кривая, приближенно отображающая характер расположения точек. По полученной кривой определяется вид приближающей функции. Рассмотрим метод наименьших квадратов, который используется для аппроксимации табличных данных с помощью непрерывной функции.
Выбираем соответствующую зависимость функции от аргумента:
y=F(x,a,b,…,c)
Неизвестными параметрами здесь являются a,b,…,c.
Приведем функцию невязок:
Для того чтобы функция невязки достигла минимального значения, необходимо потребовать равенство нулю всех ее частных производных по всем неизвестным параметрам a,b,…,c :
если полученную систему уравнений решить относительно параметров a,b,…,c, например, методом Гаусса, мы получим конкретный вид искомой функции y=F(x,a,b,…,c).
Очевидно, что изменение количества параметров не искажает сущности подхода, а лишь сказывается на количестве уравнений в системе.
Значения
найденной функции
в точках
,
,…,
будут
отличаться от значений
,
,…,
.
Значения разностей
называют отклонениями измеренных значений y от вычисленных по полученной формуле. Таким образом, можно найти величину
Среди
нескольких различных приближений одной
и той же табличной функции лучшим следует
считать то приближение, для которого
величина
принимает
наименьшее значение.
6.2 Решение задачи методом наименьших квадратов
6.2.1 Аппроксимация формулой линейной зависимости
Построим на графике заданные точки, для определения вида эмпирических формул:
Рисунок 6.2 – График функции, построенный по заданным точкам
Возьмем
в качестве аппроксимирующей формулу
линейной зависимости
.
Необходимые условия экстремума
Имеем:
П
олучаем
систему:
=>
Дальнейшие вычисления были запрограммированы на языке matlab. Текст программы приведен в приложении A.
y=0.4571*x+0.8286
Рисунок 6.3 – Результаты аппроксимации исходных данных линейной функцией (эмпирическая формула и сумма квадратов абсолютных погрешностей)
Видно,
что
достаточно велика, т.к. точки на графике
вблизи найденной прямой расположены
редко – это продемонстрировано на
следующем графике:
Рисунок 6.4 – Исходный график и график аппроксимирующей линейной функции
6.2.2 Аппроксимация формулой квадратичной зависимости
Возьмем
в качестве аппроксимирующей формулу
квадратичной зависимости
.
Необходимые условия экстремума
Нормальная система МНК:
Дальнейшие вычисления были запрограммированы на языке matlab. Текст программы приведен в приложении A.
y=-0.1476*x^2+3.4095*x-13.3429
Рисунок 6.5 – Результаты аппроксимации исходных данных квадратичной функцией (эмпирическая формула и сумма квадратов абсолютных погрешностей)
Далее приведен график исходной функции и найденной:
Рисунок 6.6 – Исходный график и график аппроксимирующей среднеквадратичной функции
Т.к.
значение
для квадратичной функции меньше значения
для линейной функции – то, полученная
эмпирическая формула
лучшим образом аппроксимирует значения исходной функции, заданной в виде набора точек
ЗАКЛЮЧЕНИЕ
В настоящее время в связи с широким распространением современных ЭВМ количественные методы исследования проникли практически во все сферы человеческой деятельности. При этом разнообразие численных методов решения как существовавших ранее, так и новых задач постоянно растет. Успех численных методов объясняется их сравнительно простой реализацией на ЭВМ. Обычное требование к численным методам состоит в том, чтобы они были устойчивыми и сходящимися. Численный метод называется устойчивым, если результаты непрерывно зависят от входных данных задачи и если погрешность округления, связанная с реализацией численного метода на ЭВМ, остается ограниченной при заданных пределах изменения параметров численного метода. Численный метод называется сходящимся, если результаты стремятся к точном решению задачи при стремлении параметров численного метода к определенным предельным значениям. Основным вопросом теории численных методов является получение численных методов, удовлетворяющих требованиям высокой точности, устойчивости и экономичности. Поэтому искусство вычислений состоит фактически не столько в предъявлении числовых результатов в виде таблиц чисел и графиков, сколько в обосновании того, что эти результаты получены с заданной точностью.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
Бахвалов Н.С. Численные методы в задачах и упражнениях. Учеб. пособие / Н.С.Бахвалов, А.В. Лапин, Е.В. Чижонков. – М.: Высш. шк., 2000.
– 190 с.
Вержбицкий В.М. Основы численных методов: Учебник для вузов / В.М. Вержбицкий. – М.: Высш. шк., 2002. – 840 с.
Доценко С.В. Конспект лекций по дисциплине «Численные методы в информатике» / С.В. Доценко. – Севастополь. Изд-во СевНТУ, 2000. – 112 с.
Приложение А
Текст программы на языке Си для численного решения дифференциальных уравнений методом Эйлера .
#include <iostream.h>
#include <conio.h>
#include <math.h>
void main() {
double a=2.0;
double b=2.5;
double h=0.05;
double x,y;
x=a;
y=1;
int i=0;
while (x<=b) { //Цикл
i++;
y=y+(2*x*pow(y,3)/(1-pow(x,2)*pow(y,3)))*h;
x+=h;
cout<<"y"<<i<<": "<<y<<endl; //Вывод у
}
cout<<endl<<"h: "<<h<<endl;
cout<<"x1: 2";
getch();
clrscr();
}
Текст программы на языке Си для численного решения дифференциальных уравнений модифицированным методом Эйлера .
#include <iostream.h>
#include <conio.h>
#include <math.h>
void main() {
float a=2;
float b=2.5;
float h=0.05;
float x,y;
x=a;
y=1;
float i=1;
cout <<"y1: "<<y<<endl;
while (x<=b) {
i++;
y=y+(2*x*pow(y,3)/(1-pow(x,2)*pow(y,3)))*h/2;
x=x+h/2;
y= y+(2*x*pow(y,3)/(1-pow(x,2)*pow(y,3)))*h;;
x=x+h/2;
cout<<"y"<<i<<": "<<y<<endl;
}
cout<<endl<<"h: "<<h<<endl;
cout<<"x1: 2";
getch();
clrscr();
}
Текст программы на языке Си для численного решения дифференциальных уравнений методом Рунге-Кутты.
#include <iostream.h>
#include <conio.h>
#include <math.h>
float f(float x, float y ); //функция
void main() {
float a=2;
float b=2.5;
float h=0.05;
float x,y,y1;
float f1,f2;
x=a;
y=1;
int i=1;
cout <<"y1: "<<y<<endl;
while (x<=b) {
i++;
f1=f(x,y);
y1=y;
y=y+f1*h;
x+=h;
f2=f(x,y);
y=y1+h*((f1+f2)/2);
cout<<"y"<<i<<": "<<y<<endl;
}
cout<<endl<<"h: "<<h<<endl;
cout<<"x1: 2";
getch();
clrscr();
}
float f(float x, float y) {
return (2*x*pow(y,3)/(1-pow(x,2)*pow(y,3)));
}
Текст программы на языке Си для нахождения экстремума, поиск отрезка методом Свенна, уточнение корня методом хорд.
#include <iostream.h>
#include <conio.h>
#include <math.h>
double x0;
double h;
double a,b; //отрезок
double x_; //искомый x*
double e=0.01;
double f(double x){
double y=15*x-2*pow(x,2)+5;
return y;
}
double f1(double x){
double y=15-4*x;
return y;
}
double f2(double x){ // вторая производая
double y=-4;
return y;
void nuton_rafson(double x0, double e){
double old_x; //предыдущее значение x
x_=x0;
do {
old_x=x_;
x_=old_x-(f1(old_x)/f2(old_x));
}
while (fabs(f1(x_))>e); // правило останова
cout<<"x*: "<<x_<<endl;
}
void svenn(int x0,int h){
int const num=10000;
double *x = new double[num];
x[0]=x0;
int k=1; //шаг 1
if (f(x[0]-h)<=f(x[0])<=f(x[0]+h)){ //шаг 2
x[1]=x[0]+h;
}
else { //шаг 3
if(f(x[0]-h)>=f(x[0])>=f(x[0]+h)){
x[1]=x[0]-h;
h=-h;
}
else {
a=x[0]-h;
b=x[0]+h;
return;
}
}
x[k+1]=x[k]+pow(2,k)*h; //шаг 4
while(f(x[k+1])>=f(x[k])){ //шаг 5
k+=1;
x[k+1]=x[k]+pow(2,k)*h;
}
if (h>0) { //шаг 6
a=x[k-1];
b=x[k+1];
}
else {
a=x[k+1];
b=x[k-1];
}
delete [] x;
}
void main(){
clrscr();
cout<<"Введите данные для метода Свенна"<<endl<<"x0=";
cin>>x0;
cout<<"h=";
cin>>h;
svenn(x0,h);
cout<<"a="<<a<<" b="<<b<<endl;
getch();
cout<<"Введите данные для метода Ньютона-Рафсона"<<endl<<"x0=";
cin>>x0;
cout<<"E=";
cin>>e;
nuton_rafson(x0,e);
cout<<"x*="<<x_;
getch();
}
Текст программы на языке Си для поиска экстремума многомерной функции методом Ньютона.
// mnogomerka.cpp: определяет точку входа для консольного приложения.
//
#include "stdafx.h"
#include <iostream>
#include <math.h>
using namespace std;
double f(double x, double y) {
return y+10*(y-sin(x))*(y-sin(x));
}
double* grad(double x, double y) {
double *a=new double[2];
a[0]=20*cos(x)*(sin(x)-y);
a[1]=-20*sin(x)+20*y+1;
return a;
}
double** gesse(double x, double y) {
double **b=new double *[2];
for (int i = 0; i < 2; i++)
b[i] = new double [2];
b[0][0]=20*(y*sin(x)+cos(2*x));
b[0][1]=-20*cos(x);
b[1][0]=-20*cos(x);
b[1][1]=-20;
return b;
}
double module (double* a) {
double s=sqrt(a[0]*a[0]+a[1]*a[1]);
return s;
}
void inversion(double **A, int N)
{
double temp;
double **E = new double *[N];
for (int i = 0; i < N; i++)
E[i] = new double [N];
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
{
E[i][j] = 0.0;
if (i == j)
E[i][j] = 1.0;
}
for (int k = 0; k < N; k++)
{
temp = A[k][k];
for (int j = 0; j < N; j++)
{
A[k][j] /= temp;
E[k][j] /= temp;
}
for (int i = k + 1; i < N; i++)
{
temp = A[i][k];
for (int j = 0; j < N; j++)
{
A[i][j] -= A[k][j] * temp;
E[i][j] -= E[k][j] * temp;
}
}
}
for (int k = N - 1; k > 0; k--)
{
for (int i = k - 1; i >= 0; i--)
{
temp = A[i][k];
for (int j = 0; j < N; j++)
{
A[i][j] -= A[k][j] * temp;
E[i][j] -= E[k][j] * temp;
}
}
}
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
A[i][j] = E[i][j];
for (int i = 0; i < N; i++)
delete [] E[i];
delete [] E;
}
double* newton(double x0,double y0, double e) {
// шаг 1
cout<<"step1"<<endl;
double *x_=new double[2],
*x_prev=new double[2],
*x=new double[2],
*d=new double[2],
t;
bool prev_yes=false; //выполнилось ли при x_(k-1)
x_prev[0]=x[0]=x0;
x_prev[1]=x[1]=y0;
double fx=f(x[0],x[1]);
do{
// шаг 3
cout<<"step3"<<endl;
// Находим градиент функции
double* g = grad(x[0],x[1]);
// шаг 4
cout<<"step4"<<endl;
if(module(g)<=e){
x_[0]=x[0];
x_[1]=x[1];
return x_;
}
else{ // шаг 6
cout<<"step6"<<endl;
// находим матрицу Гессе
double **h=new double *[2];
for (int i = 0; i < 2; i++)
h[i] = new double [2];
h=gesse(x[0],x[1]);
// шаг 7
cout<<"step7"<<endl;
// находим обратную матрицу
inversion(h,2);
// шаг 8
cout<<"step8"<<endl;
if(h[0][0]>0&&(h[0][0]*h[1][1]-h[0][1]*h[1][0])>0){
// шаг 9
cout<<"step9"<<endl;
d[0]=-h[0][0]*g[0]+(-h[0][1])*g[1];
d[1]=-h[1][0]*g[0]+(-h[1][1])*g[1];
t=1;
}
else{
d[0]=-g[0];
d[1]=-g[1];
if(f(x[0],x[1])<f(x_prev[0],x_prev[1])) t=-1;
else t=1;
}
// шаг 10
cout<<"step10"<<endl;
x_prev[0]=x[0];
x_prev[1]=x[1];
x[0]=x[0]+t*d[0];
x[1]=x[1]+t*d[1];
//шаг 11
cout<<"step11"<<endl;
d[0]=x[0]-x_prev[0];
d[1]=x[1]-x_prev[1];
if(module(d)<e&&fabs(f(x[0],x[1])-f(x_prev[0],x_prev[1]))<e&&prev_yes){
x_[0]=x[0];
x_[1]=x[1];
return x_;
}
else if(module(d)<e&&fabs(f(x[0],x[1])-f(x_prev[0],x_prev[1]))<e)
prev_yes=true;
}
}
while(1);
}
int _tmain(int argc, _TCHAR* argv[])
{
double *x=new double[2];
x=newton(1,1,0.1);
cout<<"Минимум функции достигается при"<<endl;
cout<<"x = "<<x[0]<<endl;
cout<<"y = "<<x[1]<<endl;
cout<<"и имеет значение f = "<<f(x[0],x[1])<<endl;
delete [] x;
std::cin.get();
std::cin.get();
return 0;
}
Текст программы на языке matlab для метода наименьших квадратов (методы обработки экспериментальных данных).
function y = LinearFunc(k, x, b)
%Функция y = kx + b
y = k*x + b;
function K = LinearApproximation(X, Y, size)
%---------------------------Аппроксимация формулой линейной зависимости
fprintf(' Аппроксимация линейной зависимостью\n');
A = zeros(2, 2); %Коэффициенты при неизвестных (матрица)
B = zeros(2, 1); %Вектор-столбец правых частей
%Заполнение матрицы коээфициентов и столбца правых частей
for i = 1 : size
A(1, 1) = A(1, 1) + X(i)*X(i);
A(1, 2) = A(1, 2) + X(i);
B(1, 1) = B(1, 1) + X(i)*Y(i);
B(2, 1) = B(2, 1) + Y(i);
end
A(2, 1) = A(1, 2);
A(2, 2) = size;
%Вывод матрицы для проверки (если необходимо)
fprintf(' Расширенная матрица системы уравнений:\n');
for i = 1 : 2
fprintf(' ');
for j = 1 : 2
fprintf('%.6f ', A(i, j));
end
fprintf('%.6f\n', B(i, 1));
end
%Решение системы уравнений методом обратной матрицы
K = inv(A)*B; %Вектор, содержащий решение
End
function y = QuadFunc(a, b, c, x)
%Вычисление функции y = ax^2 + bx + c
y = a*(x^2) + b*x + c;
function K = QuadraticApproximation(X, Y, size)
%-------------------------Аппроксимация формулой квадратичной зависимости
fprintf(' Аппроксимация квадратичной зависимостью\n');
A = zeros(3, 3); %Коэффициенты при неизвестных (матрица)
B = zeros(3, 1); %Вектор-столбец правых частей
%Заполнение матрицы коээфициентов и столбца правых частей
for i = 1 : size
A(1, 1) = A(1, 1) + X(i)^4;
A(1, 2) = A(1, 2) + X(i)^3;
A(1, 3) = A(1, 3) + X(i)^2;
A(2, 3) = A(2, 3) + X(i);
B(1, 1) = B(1, 1) + Y(i)*(X(i)^2);
B(2, 1) = B(2, 1) + X(i)*Y(i);
B(3, 1) = B(3, 1) + Y(i);
end
A(2, 1) = A(1, 2);
A(2, 2) = A(1, 3);
A(3, 1) = A(2, 2);
A(3, 2) = A(2, 3);
A(3, 3) = size;
%Вывод матрицы для проверки (если необходимо)
fprintf(' Расширенная матрица системы уравнений:\n');
for i = 1 : 3
fprintf(' ');
for j = 1 : 3
fprintf('%.6f ', A(i, j));
end
fprintf('%.6f\n', B(i, 1));
end
%Решение системы уравнений методом обратной матрицы
K = inv(A)*B; %Вектор-решение
end