p align="left">Рішення. Умову завдання перепишемо наступним чином . Приймаючи що і то коротко систему рівнянь можна записати так . Якщо відомо деяке наближення кореня системи рівнянь, то поправку можна знайти рішаючи систему . Розкладемо ліву частину рівняння по степеням малого вектору , обмежуючись лінійними членами . == - матриця похідних (матриця Якобі) (). Складемо матрицю похідних (матрицю Якобі): Якщо , то , де - матриця обернена до матриці Якобі. Таким чином послідовне наближення кореня можна обчислити за формулою або . Умовою закінчення ітераційного процесу наближення корення вибираємо умову , - евклідова відстань між двома послідовними наближеннями ;- число, що задає мінімальне наближення. Для рішення систем нелінійних рівнянь за даним алгоритмом призначена програма Work3.cpp //------------------------------------------------------------ // Work3.cpp //------------------------------------------------------------ // "Числові методи" // Завдання 3 // Розв'язування системи нелінійних рівнянь методом Ньютона #include <stdio.h> #include <iostream.h> #include <conio.h> #include <math.h> #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 <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]. Рішення. В загальному вигляді ряд Фурьє функції виглядає так: , де =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
Страницы: 1, 2, 3
|