заданим значенням X[] елементи матриці W ;
Вихідні дані:
X[N] - вектор наближеного значення мінімуму.
Функція повертає код помилки
0 - система рівнянь успішно розв’язана
1 - det W=0 */
{int n=N;
float len;
float W[N][N],Winv[N][N],Y[N],deltaX[N];
do
{pJakobi(X,W);
if(invMatr(n,W,Winv)) return 1;
pSist(X,Y);
DobMatr(n,Winv,Y,deltaX);
X[0]-=deltaX[0];
X[1]-=deltaX[1];
len=sqrt(deltaX[0]*deltaX[0]+deltaX[1]*deltaX[1]);}
while (len>eps);
return 0;}
//int main()
{float X[N],eps;
// початкові умови
eps=.0001;
X[0]=0.0; X[1]=1.0;
if (NelinSist(X,fJakobi,fSist,eps))
{ cout << "Error of matrix: detW=0"; return 1;}
printf("X= %5.4f Y= %5.4f\n",X[0],X[1]);
cout << "\n Press any key ...";
getch();}
Результат роботи програми:
X= 0.1477 Y= 1.0214
Завдання 4
Знайти точку мінімуму та мінімальне значення функції
методом Ньютона.
Рішення.
Матриця Гессе
Ітераційний процес послідовного наближення мінімуму функції буде таким:
де
Для закінчення ітераційного процесу використаємо умову
Для пошуку мінімуму функції за методом Ньютона призначена програма Work4.cpp
//------------------------------------------------------------
// matrix.h
//-----------------------------------------------------------
const int nMax=2; // кількість рівнянь
const float ZERO=.00000001;
int invMatr(int n,float A[nMax][nMax],float Ainv[nMax][nMax])
/* Функція знаходить обернену матрицю
Вхідні дані:
A- масив з коефіцієнтами при невідомих;
n- порядок матриці А(кількість рівнянь системи);
Вихідні дані:
Ainv- матриця обернена до матриці А;
функція повертає код помилки:
0- помилки немає;
1- матриця А вироджена. */
{float aMax,t; // максимальний елемент , тимчасова змінна
int i,j,k,l;
// формуємо одиничну матрицю
for(i=0; i<n; i++)
for (j=0; j<n; j++)
Ainv[i][j] = (i==j)? 1. : 0.;
for (k=0; k<n; k++)
{// знаходимо мах по модулю елемент
aMax=A[k][k]; l=k;
for (i=k+1; i<n; i++)
if (fabs(A[i][k])>fabs(aMax))
{ aMax=A[i][k]; l=i; }
// якщо модуль головного елементу aMax менший за програмний 0 (ZERO)
if ( fabs(aMax)<ZERO ) return 1;
// якщо потрібно, міняємо місцями рівняння Pk i Pl
if ( l!=k)
for( j=0; j<n; j++)
{t=A[l][j]; A[l][j]=A[k][j]; A[k][j]=t;
t=Ainv[l][j]; Ainv[l][j]=Ainv[k][j]; Ainv[k][j]=t;}
// ділимо k-й рядок на головний елемент
for (j=0; j<n; j++) { A[k][j]/=aMax; Ainv[k][j]/=aMax; }
// обчислюємо елементи решти рядків
for (i=0; i<n; i++)
if( i!=k )
{t=A[i][k];
for (j=0; j<n; j++)
{A[i][j]-=t*A[k][j];
Ainv[i][j]-=t*Ainv[k][j];}}}
return 0;}
void DobMatr(int n, float A[nMax][nMax], float B[nMax],float X[nMax])
// функція знаходить добуток матриці А на вектор В і результат повертає в
// векторі Х
{int i,j;
float summa;
for (i=0; i<n; i++)
{summa=0;
for (j=0; j<n; j++)
{summa+=A[i][j]*B[j];
X[i]=summa;}}
} // DobMatr
//------------------------------------------------------------
// Work4.cpp
//------------------------------------------------------------
// "Числові методи"
// Завдання 4
// Пошук мінімуму функції методом Ньютона
#include <stdio.h>
#include <iostream.h>
#include <conio.h>
#include <math.h>
#include "matrix.h"
const int N=2; // степінь матриці Гессе
float myFunc(float x[N])
{ return exp(-x[1])-pow(x[1]+x[0]*x[0],2); }
typedef void (*funcH) (float[N], float[N][N]);
void fHesse(float X[N],float H[N][N])
// функції , які складають матрицю Гессе
{H[0][0]=-4.*X[1]-6.*X[0]*X[0]; H[0][1]=-4.*X[0];
H[1][0]=-4; H[1][1]=exp(-X[1])-21;}
typedef void (*funcG) (float[N], float[N]);
void fGrad(float X[N],float Y[N])
{Y[0]=-4*X[1]*X[0]-3*X[0]*X[0]*X[0];
Y[1]=exp(-X[1])-2.*X[1]-2*X[0]*X[0];}
//int fMin(float X[N], funcG pGrad, funcH pHesse,float eps)
/* Функція знаходить точку мінімуму рівняння методом Ньютона.
Вхідні дані:
X[N] - вектор значень початкового наближення
pGrad - вказівник на функцію, яка обчислює по
заданим значенням X[] значення grad f(X) ;
pHesse - вказівник на функцію, яка обчислює по
заданим значенням X[] елементи матриці H ;
Вихідні дані:
X[N] - вектор наближеного значення мінімуму.
Функція повертає код помилки
0 - система рівнянь успішно розв’язана
1 - det H=0 */
{int n=N;
float modGrad;
float Hesse[N][N],HesseInv[N][N],Grad[N],deltaX[N];
do
{pHesse(X,Hesse);
if(invMatr(n,Hesse,HesseInv)) return 1;
pGrad(X,Grad);
DobMatr(n,HesseInv,Grad,deltaX);
X[0]-=deltaX[0];
X[1]-=deltaX[1];
modGrad=sqrt(deltaX[0]*deltaX[0]+deltaX[1]*deltaX[1]);}
while (modGrad>eps);
return 0;}
//int main()
{float X[N],eps;
// початкові умови
eps=.0001;
X[0]=0.5; X[1]=0.5;
if (fMin(X,fGrad,fHesse,eps))
{ cout << "Error of matrix: detH=0"; return 1;}
printf("X= %5.5f Y= %5.4f\n f(x,y)= %4.3f\n ",X[0],X[1],myFunc(X));
cout << "\n Press any key ...";
getch();}
Результат роботи програми:
x= -0.0000 y= 0.3523
f(x,y)= 0.579
Завдання 5
Розкласти в ряд Фурьє функцію
Рішення.
В загальному вигляді ряд Фурьє функції
В нашому випадку відрізок розкладення функції – [-1; 1], тому проводимо лінійну заміну змінної
Для наближеного обчислення коефіцієнтів ряду Фурьє використаємо квадратурні формули, які утворюються при інтерполяції алгебраїчним многочленом підінтегральних функцій
де
Для обчислення наближених значень коефіцієнтів ряду Фурьє по формулам (1), (2), (3) призначена процедура (функція) Fourier.
//---------------------------------------------------------
// Work5.h
//---------------------------------------------------------
#include <math.h>
const double Pi=3.141592653;
// функція повертає і-й вузол квадратурної формули, 2N+1-кілікість вузлів
inline double FuncXi(int N, int i) {return -Pi+(2*Pi*i)/(2*N+1);}
typedef double (*Func)(double); // опис типу вказівника на функцію
char Fourier(Func F_name, int CountN, int CountK,double **Arr)
/* функція обчислює коефіцієнти ряду Фурьє
Вхідні дані:
F_mame - вказівник на функцію(функтор), яка обчислює значення функції
f(x) на відрізку [-п; п];
CountN - число, яке задає розбиття відрізка [-п; п] на рівні частини
довжиною 2п/(2*CountN+1);
CountK - кількість обчислюваних пар коефіцієнтів;
Вихідні дані:
Arr - двомірний масив розміру [CountK+1][2], в якому
знаходяться обчислені коефіцієнти ряду Фурьє.
Функція повертає значення коду помилки:
Fourier=0 - помилки немає;
Fourier=1 - якщо CountN<CountK;
Fourier=2 - якщо CountK<0;*/
{double a,b,sumA,sumB;
int i,k;
if (CountN < CountK) return 1;
if (CountK < 0) return 2;
// обчислення а0
sumA=0;
for (i=0; i< 2*CountN+1; i++) sumA+=F_name(FuncXi(CountN,i));
a=1./(2*CountN+1)*sumA;
Arr[0][0]=a;
// обчислення коефіцієнтів аk,bk
for (k=1; k<=CountK; k++)
{sumA=sumB=0;
for (i=0; i<2*CountN+1; i++)
{sumA+=F_name(FuncXi(CountN,i))*cos(2*Pi*k*i/(2*CountN+1));
sumB+=F_name(FuncXi(CountN,i))*sin(2*Pi*k*i/(2*CountN+1));}
a=(2./(2*CountN+1))*sumA;
b=(2./(2*CountN+1))*sumB;
Arr[k][0]=a;
Arr[k][1]=b;}
return 0;}
//------------------------------------------------------------
// Work5.cpp
//------------------------------------------------------------
// "Числовы методи"
// Завдання 5
// Розрахунок коэфіцієнтів ряду Фурьє
#include "Work5.h"
#include <stdio.h>
#include <iostream.h>
#include <conio.h>
double f(double x)
// функція повертає значення функції f(x)
{return sqrt(Pi*Pi*x*x+1);}
const int N=20; // константа, яка визначає розбиття відрізка [-п; п]
// на рівні частини
const int CountF=15; // кількість пар коефіцієнтів ряду
void main()
{double **data;
data = new double *[CountF+1];
for ( int i=0; i<=CountF; i++) data[i] = new double [2];
if (Fourier(f,N,CountF,data) != 0)
{cout << "\n Помилка !!!";
return;}
// Вивід результатів
printf("a0= %lf\n",data[0][0]);
for (int i=1;i<=CountF;i++)
printf("a%d = %lf , b%d = %lf\n",i,data[i][0],i,data[i][1]);
cout << " Press any key ...";
getch();}
Результат роботи програми Work5.cpp
|