Числові методи (стр. 4 из 4)

заданим значенням 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];

{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];

{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].

Рішення.

В загальному вигляді ряд Фурьє функції

виглядає так:

, де

=0, 1, 2, …

В нашому випадку відрізок розкладення функції – [-1; 1], тому проводимо лінійну заміну змінної

. Тоді умова завдання стане такою:

Для наближеного обчислення коефіцієнтів ряду Фурьє використаємо квадратурні формули, які утворюються при інтерполяції алгебраїчним многочленом підінтегральних функцій

(1)

(2)

(3)

де

– число вузлів квадратурної формули;

– вузли квадратурної формули ,

=0, 1, 2, …, 2

Для обчислення наближених значень коефіцієнтів ряду Фурьє по формулам (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