Числові методи
Контрольная работа - Компьютеры, программирование
Другие контрольные работы по предмету Компьютеры, программирование
перепишемо наступним чином
.
Приймаючи що і то коротко систему рівнянь можна записати так
.
Якщо відомо деяке наближення кореня системи рівнянь, то поправку можна знайти рішаючи систему
.
Розкладемо ліву частину рівняння по степеням малого вектору , обмежуючись лінійними членами
.
== матриця похідних (матриця Якобі) ().
Складемо матрицю похідних (матрицю Якобі):
Якщо , то ,
де матриця обернена до матриці Якобі.
Таким чином послідовне наближення кореня можна обчислити за формулою
або
.
Умовою закінчення ітераційного процесу наближення корення вибираємо умову
,
евклідова відстань між двома послідовними наближеннями ; число, що задає мінімальне наближення.
Для рішення систем нелінійних рівнянь за даним алгоритмом призначена програма
Work3.cpp
//------------------------------------------------------------
// Work3.cpp
//------------------------------------------------------------
// "Числові методи"
// Завдання 3
// Розвязування системи нелінійних рівнянь методом Ньютона
#include
#include
#include
#include
#include "matrix.h"
const int N=2; // степінь матриці Якобі (кількість рівнянь)
typedef void (*funcJ) (float[N], float[N][N]);
void fJakobi(float X[N],float J[N][N])
// функції , які складають матрицю Гессе
{J[0][0]=cos(X[0]); J[0][1]=cos(X[1]);
J[1][0]=2*X[0]; J[1][1]=-2*X[1]+1;}
typedef void (*funcF) (float[N], float[N]);
void fSist(float X[N],float Y[N])
{Y[0]=sin(X[0])+sin(X[1])-1;
Y[1]=X[0]*X[0]-X[1]*X[1]+X[1];}
//int NelinSist(float X[N], funcJ pJakobi, funcF pSist,float eps)
/* Функція знаходить кореня системи нелінійних рівнянь методом Ньютона.
Вхідні дані:
X[N] - вектор значень початкового наближення
pSist - вказівник на функцію, яка обчислює по
заданим значенням X[] значення функції f(X) ;
pJakobi - вказівник на функцію, яка обчислює по
заданим значенням 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
#include
#include
#include
#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[] елемен?/p>