Числові методи

МІНІСТЕРСТВО ОСВІТИ УКРАЇНИ


ЧЕРНІВЕЦЬКИЙ ДЕРЖАВНИЙ УНІВЕРСИТЕТ


ІМ. Ю. ФЕДЬКОВИЧА


КОНТРОЛЬНА РОБОТА


з дисципліни " Числові методи "


Варіант 16.


Виконав


студент 2-го курсу


кафедри ЕОМ


Перевірив


м. Чернівці


Завдання 1


Задана СЛАР



а) розв’язати цю систему методом Гауса за схемою з частковим вибором головного елементу;


б)розв’язати цю систему за формулою


.


– вектор невідомих, – вектор вільних членів, – обернена матриця до матриці з коєфіцієнтів при невідомих.


Обернену матрицю знай ти методом Гауса - Жордана за схемою з частковим вибором головного елемента.


Рішення.


а) Прямий хід методу Гауса.


()


Запишемо матрицю .



1-й крок.


Серед елементів першого стовпчика шукаємо максимальний:



Перше і друге рівняння міняємо місцями.



Розділимо рівняння (1) на 2.5


(1)


Від рівняння (2) віднімемо 1.7Р1 .




(2)




(3)


Таким чином в кінці першого кроку отримуємо систему



2-й крок.


Порядок рівнянь зберігається.


(2)




(3)


Після другого кроку система рівнянь стала такою:



Зворотній хід.


З рівняння (3) ;


з рівняння (2) ;


з рівняння (1) ;


Для рішення системи лінійних рівнянь методом Гауса призначена програма Work1_1.


//------------------------------------------------------------


// Work1_1.cpp


//------------------------------------------------------------


// "Числові методи"


// Завдання 1


// Рішення системи лінійних рівнянь методом Гауса


#include <stdio.h>


#include <iostream.h>


#include <conio.h>


const int nMax=5; // максимальна кількість рівнянь


const float ZERO=.0000001;


int fGaus(float A[nMax][nMax],float B[nMax],int n,float X[nMax])


/* Функція розв'язує систему лінійних рівнянь методом Гауса за схемою з


частковим вибором головного елементу.


Вхідні дані:


A- масив з коефіцієнтами при невідомих;


В- масив з вільними членами СЛАР;


n- порядок матриці А(кількість рівнянь системи);


Вихідні дані:


Х- масив з коренями системи;


функція повертає код помилки:


0- сисетма успішно розв’язана;


1- матриця А вироджена. */


{float aMax,t; // максимальний елемент , тимчасова змінна


int i,j,k,l;


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=B[l]; B[l]=B[k]; B[k]=t;}


// ділимо k-те рівняння на головний елемент


for (j=0; j<n; j++) A[k][j]/=aMax;


B[k]/=aMax;


// обчислюємо коефіцієнти A[i][j] та вільні члени решти рівнянь


for (i=k+1; i<n; i++)


{t=A[i][k]; B[i]-=t*B[k];


for (j=0; j<n; j++) A[i][j]-=t*A[k][j];}


} // for (k)


// Зворотній хід


for ( k=n-1; k>=0; k--)


{X[k]=0;


for (l=k+1; l<n; l++) X[k]+=A[k][l]*X[l];


X[k]=B[k]-X[k];}


return 0;


} // fGaus()


void main()


{float A[nMax][nMax];


float B[nMax];


float X[nMax];


int n,i,j;


char *strError="n Error of file !";


FILE *FileIn,*FileOut;


FileIn=fopen("data_in.txt","r"); // відкриваємо файл для читання


if (FileIn==NULL)


{cout << " "Data_in.txt": Error open file or file not found !!!n";


goto exit;}


FileOut=fopen("data_out.txt","w"); // відкриваємо файл для запису


if (FileOut==NULL)


{cout << " "Data_out.txt": Error open file !!!n";


goto exit;}


if(fscanf(FileIn,"%d",&n)==NULL)


{ cout << strError; goto exit;};


for (i=0; i<n; i++)


for(j=0; j<n; j++)


fscanf(FileIn,"%f",&(A[i][j]));


for (i=0; i<n;i++)


if(fscanf(FileIn,"%f",&(B[i]))==NULL)


{ cout << strError; goto exit;}


if(fGaus(A,B,n,X)!=0)


{ cout << "n det|A|=0 !"; goto exit;}


// Вивід результатів


for (i=0; i<n; i++)


{printf(" x[%d]= %f ",i+1,X[i]);


fprintf(FileOut," x[%d]= %f ",i+1,X[i]);}


fclose(FileIn);


fclose(FileOut);


exit: cout << "n Press any key ...";


getch();}


Результат роботи програми:


x[1]= 3.017808 x[2]= 0.356946 x[3]= -0.302131


б) Знайдемо обернену матрицю .


0-й крок.


А Е


1-й крок.






;



2-й крок.




;



3-й крок.


; ;






.


Даний алгоритм рішення системи лінійних рівнянь реалізований в програмі Work1_2.


//------------------------------------------------------------


// Work1_2.cpp


//------------------------------------------------------------


// "Числові методи"


// Завдання 1


// Рішення системи лінійних рівнянь методом Гауса-Жордана


#include <stdio.h>


#include <iostream.h>


#include <conio.h>


const int nMax=5; // максимальна кількість рівнянь


const float ZERO=.0000001;


int fGausJordan(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;


} // fGausJordana()


void fDobMatr(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;}}


} // fDobMatr


void main()


{float A[nMax][nMax],Ainv[nMax][nMax];


float B[nMax];


float X[nMax];


int n,i,j;


char *strError="n Error of file !";


FILE *FileIn,*FileOut;


FileIn=fopen("data_in.txt","r"); // відкриваємо файл для читання


if (FileIn==NULL)


{cout << " "Data_in.txt": Error open file or file not found !!!n";


goto exit;}


FileOut=fopen("data_out.txt","w"); // відкриваємо файл для запису


if (FileOut==NULL)


{cout << " "Data_out.txt": Error open file !!!n";


goto exit;}


if(fscanf(FileIn,"%d",&n)==NULL)


{ cout << strError; goto exit;};


for (i=0; i<n; i++)


for(j=0; j<n; j++)


fscanf(FileIn,"%f",&(A[i][j]));


for (i=0; i<n;i++)


if(fscanf(FileIn,"%f",&(B[i]))==NULL)


{ cout << strError; goto exit;}


if(fGausJordan(n,A,Ainv)!=0)


{ cout << "n det|A|=0 !"; goto exit;}


fDobMatr(n,Ainv,B,X);


// Вивід результатів


for (i=0; i<n; i++)


{printf(" x[%d]= %f ",i+1,X[i]);


fprintf(FileOut," x[%d]= %f ",i+1,X[i]);}


fclose(FileIn);


fclose(FileOut);


exit: cout << "n Press any key ...";


getch();}


Результат роботи програми:


x[1]= 3.017808 x[2]= 0.356946 x[3]= -0.302131


Завдання 2


Задана задача Коші


,


а) Знайти розв’язок в табличній формі методом Рунге-Кутта:

















, , .


б) Інтерполювати цю функцію кубічним сплайном. Систему рівнянь для моментів кубічного сплайну розв’язати методом прогонки. Вибрати крайові умови для кубічного сплайну у вигляді


.


в) Використовуючи кубічний сплайн з пункту б) обчислити методом Сімпсона .


Взяти (– кількість відрізків розбиття).


Рішення.


а) Метод Рунге-Кутта


Розрахунок будемо проводити за наступними формулами :


;


;


;


;


;


.


Цей алгоритм реалізовується в програмі Work2_1.


//------------------------------------------------------------


// Work2_1.cpp


//------------------------------------------------------------


// "Числові методи"


// Завдання 2


// Рішення задачі Коші методом Рунге-Кутта


#include <stdio.h>


#include <iostream.h>


#include <conio.h>


typedef float (*pfunc)(float,float); // pfunc - вказівник на функцію


const int nMax=5; // максимальна кількість відрізків розбиття


void fRunge_Kutta(pfunc f, float x0, float y0,float h, int n, float Y[nMax])


/* Функція знаходить табличне значення функції методом Рунге-Кутта


Вхідні дані:


f - функція f(x,y)


x0,y0 - початкова точка;


h - крок;


n- кількість точок розбиття;


Вихідні дані:


Y- вектор значень функції*/


{float k1,k2,k3,k4,x; // максимальний елемент , тимчасова змінна


int i;


x=x0; Y[0]=y0;


for (i=0; i<n-1; i++)


{k1=f(x,Y[i]);


k2=f(x+h/2, Y[i]+k1*h/2);


k3=f(x+h/2, Y[i]+k2*h/2);


k4=f(x+h, Y[i]+h*k3);


Y[i+1]=Y[i]+(h/6)*(k1+2*k2+2*k3+k4);


x+=h;}}


float Myfunc(float x,float y)


{return log10(cos(x+y)*cos(x+y)+2)/log10(5);}


void main()


{float Y[nMax],h,x0,y0;


int n,i;


char *strError="n Error of file !";


FILE *FileIn,*FileOut, *FileOut2;


FileIn=fopen("data2_in.txt","r"); // відкриваємо файл для читання


if (FileIn==NULL)


{cout << " "Data2_in.txt": Error open file or file not found !!!n";


goto exit;}


FileOut=fopen("data2_out.txt","w"); // відкриваємо файл для запису


if (FileOut==NULL)


{cout << " "Data2_out.txt": Error open file !!!n";


goto exit;}


FileOut2=fopen("data2_2in.txt","w"); // відкриваємо файл для запису


if (FileOut==NULL)


{cout << " "Data2_2in.txt": Error open file !!!n";


goto exit;}


if(fscanf(FileIn,"%d%f%f%f,",&n,&h,&x0,&y0)==NULL)


{ cout << strError; goto exit;};


fRunge_Kutta(Myfunc,x0,y0,h,n,Y);


// Вивід результатів


for (i=0; i<n; i++)


{printf(" x[%d]= %4.2f ",i,Y[i]);


fprintf(FileOut," x[%d]= %4.2f ",i,Y[i]);


fprintf(FileOut2,"%4.2f ",Y[i]);}


fclose(FileIn);


fclose(FileOut);


exit: cout << "n Press any key ...";


getch();}


Результат роботи програми (файл "data2_out.txt"):


x[0]= 1.00 x[1]= 1.05 x[2]= 1.10 x[3]= 1.14 x[4]= 1.18


б) В загальному вигляді кубічний сплайн виглядає наступним чином:


,


Параметри кубічного сплайну будемо обчислювати , використовуючи формули:


; ;


; , де


– моменти кубічного сплайну.


Моменти мають задовольняти такій системі рівнянь:


.


Для ; ; .


Якщо прийняти до уваги граничні умови , то систему можна записати так


.


В даному випадку матриця з коефіцієнтів при невідомих є тридіагональною


,


тому для знаходження моментів кубічних сплайнів застосуємо метод прогонки.


На прямому ході обчислюємо такі коефіцієнти.


; ;


На зворотньому ході обчислюємо значення моментів кубічного сплайну.


; .


Для знаходження коефіцієнті вкубічного сплайну призначена програма Work2_2.


//------------------------------------------------------------


// Work2_2.cpp


//------------------------------------------------------------


// "Числові методи"


// Завдання 2


// Інтерполювання функції кубічним сплайном


#include <stdio.h>


#include <iostream.h>


#include <conio.h>


const int nMax=4; // максимальна кількість відрізків розбиття


const float x0=0.;// початкова точка сітки


const float h=0.1;// крок розбиття


// вектори матриці А


float a[]={0., 0.5, 0.5};


float b[]={2., 2., 2.};


float c[]={0.5, 0.5, 0.};


//void fMetodProgonku( int n,float a[nMax],float b[nMax],float c[nMax],float d[nMax], float M[nMax+1])


/* Функція знаходить моменти кубічного сплайну методом прогонки


Вхідні дані:


a,b,c -вектори матриці А ;


d - вектор вільних членів;


n- степінь матриці А;


Вихідні дані:


М- вектор моментів кубічного сплайну.*/


{float k[nMax],fi[nMax];


int i;


// прямий хід


for (i=0; i<n; i++)


{k[i] = (i==0)? -c[i]/b[i] : -c[i]/(b[i]+a[i]*k[i-1]);


fi[i] = (i==0)? d[i]/b[i] : (-a[i]*fi[i-1]+d[i])/(b[i]+a[i]*k[i-1]);}


//зворотній хід


for (i=n; i>0; i--)


M[i] = (i==n)? fi[i-1] : k[i-1]*M[i+1]+fi[i-1];}


void fSplain( int n,float h,float Y[nMax+1],float M[nMax+1],float Ak[nMax][4])


/* Функція обчислює коефіцієнти кубічного сплайну


Вхідні дані:


n- кількість відрізків розбиття;


H - крок розбиття відрізку [X0; Xn]]


М- вектор моментів кубічного сплайну.


Y- вектор значень функції f(x,y) в точках x[0],x[1],...x[n].


Вихідні дані:


Ak- матриця коефіцієнтів кубічного сплайну.*/


{int i;


for (i=0; i<n; i++)


{Ak[i][0] = Y[i];


Ak[i][1] = (Y[i+1]-Y[i])/h-h/6*(2.*M[i]+M[i+1]);


Ak[i][2] = M[i]/2;


Ak[i][3] = (M[i+1]-M[i])/6*h;}}


void main()


{float Y[nMax+1],d[nMax],M[nMax+1],Ak[nMax][4];


int n,i,j;


n=nMax;


M[0]=0; M[n]=0; //крайові умови


char *strError="n Error of file !";


FILE *FileIn,*FileOut,*FileOut2;


FileIn=fopen("data2_2in.txt","r"); // відкриваємо файл для читання


if (FileIn==NULL)


{ cout << " "Data2_2in.txt": Error open file or file not found !!!n";


goto exit; }


FileOut=fopen("data2_2ou.txt","w"); // відкриваємо файл для запису


if (FileOut==NULL)


{ cout << " "Data2_2ou.txt": Error open file !!!n";


goto exit; }


FileOut2=fopen("data2_3in.txt","w"); // відкриваємо файл для запису


if (FileOut2==NULL)


{ cout << " "Data2_3in.txt": Error open file !!!n";


goto exit; }


// читаємо вектор Y


for (i=0; i<=n; i++)


if(fscanf(FileIn,"%f,",&(Y[i]))==NULL)


{ cout << strError; goto exit;};


// обчислюємо вектор d


for (i=1; i<n; i++) d[i-1]=3/(h*h)*(Y[i+1]-2*Y[i]+Y[i-1]);


//fMetodProgonku(n-1,a,b,c,d,M);


fSplain( n,h,Y,M,Ak);


// Вивід результатів в тому числі і для наступного завдання


fprintf(FileOut2,"%dn",n); // n - кількість відрізків


// координати точок сітки по Х


for(float xi=x0,i=0; i<n; i++) fprintf(

FileOut2,"%2.2f ",xi+h*i);


fprintf(FileOut2,"n");


for (i=0; i<n; i++)


{for (j=0; j<4; j++)


{printf("a[%d,%d]= %4.4f ",i,j,Ak[i][j]);


fprintf(FileOut,"a[%d,%d]= %4.4f ",i,j,Ak[i][j]);


fprintf(FileOut2,"%4.4f ",Ak[i][j]);}


cout << endl;


fprintf(FileOut,"n");


fprintf(FileOut2,"n");}


fclose(FileIn);


fclose(FileOut);


exit: cout << "n Press any key ...";


getch();}


Результат роботи програми (" data2_2uo.txt"):


a[0,0]= 1.0000 a[0,1]= 0.5104 a[0,2]= 0.0000 a[0,3]= -0.0104


a[1,0]= 1.0500 a[1,1]= 0.4793 a[1,2]= -0.3107 a[1,3]= 0.0118


a[2,0]= 1.0960 a[2,1]= 0.4525 a[2,2]= 0.0429 a[2,3]= -0.0068


a[3,0]= 1.1410 a[3,1]= 0.4407 a[3,2]= -0.1607 a[3,3]= 0.0054


в) Розіб’ємо відрізок на частин.


Складова формула Сімпсона буде мати вигляд:


;


де - крок розбиття, – значення функції в точках сітки.


Замінимо значеннями кубічних сплайнів із пункту б) цього завдання.



Для оцінки похибки використаємо правило Рунге. Для цього обчислимо наближені значення інтегралу з кроком (), а потім з кроком ().


За наближене значення інтегралу, обчисленого за формулою Сімпсона з поправкою по Рунге приймемо: .


//------------------------------------------------------------


// Work2_3.cpp


//------------------------------------------------------------


// "Числові методи"


// Завдання 2


// Обчислення інтегралу методом Сімпсона з використанням кубічного сплайну


#include <stdio.h>


#include <iostream.h>


#include <conio.h>


#include <math.h>


// визначення сплайнового класу


class Tsplain


{public:


int kol; // кількість рівнянь (відрізків розбиття)


float ** Ak; // масив коефіцієнтів


float * Xi; // вектор початків відрізків


float vol(float x); // функція повертає значення сплайну в точці х


Tsplain(int k); // constructor};


Tsplain::Tsplain(int k)


{kol=k;


Xi=new float[kol];


Ak=new float*[kol];


for(int i=0; i<kol; i++) Ak[i]=new float[kol];}


float Tsplain::vol(float x)


{float s=0.;


int i,t;


// шукаємо відрізок t де знаходиться точка х


for (i=0; i<kol; i++) if (x>=Xi[i]) { t=i; break; }


s=Ak[t][0];


for (i=1; i<kol; i++)


s+=Ak[t][i]*pow((x-Xi[t]),i);


return s;}


float fSimps(float down,float up, int n, Tsplain *spl)


/* Функція обчислює інтеграл методом Сімпсона з використанням кубічного сплайну


Вхідні дані:


down,up -границі інтегрування ;


n- число відрізків , на яке розбиваєтьться відрізок інтегрування ;


spl - вказівник на об’єкт класу Tsplain ( кубічний сплайн );


Вихідні дані:


функція повертає знайдене значення інтегралу.*/


{float s=0;


float h=(up-down)/(float)n;


int i;


s=spl->vol(down)+spl->vol(up-h);


for (i=2; i<n; i+=2)


s+=2*(spl->vol(down+i*h));


for (i=1; i<n; i+=2)


s+=4*(spl->vol(down+i*h));


return s*h;}


void main()


{int kol; // кількість рівняннь кубічного сплайну


float down,up;


float I1,I2,I,eps;


int n,i,j;


char *strError="n Error of file !";


FILE *FileIn,*FileOut;


FileIn=fopen("data2_3in.txt","r"); // відкриваємо файл для читання


if (FileIn==NULL)


{ cout << " "Data2_3in.txt": Error open file or file not found !!!n";


goto exit; }


FileOut=fopen("data2_3ou.txt","w"); // відкриваємо файл для запису


if (FileOut==NULL)


{ cout << " "Data2_3ou.txt": Error open file !!!n";


goto exit; }


// читаємо kol


if(fscanf(FileIn,"%d,",&kol)==NULL)


{ cout << strError; goto exit;};


Tsplain *sp;


sp=new Tsplain(kol);


// читаємо вектор Xi


for(i=0; i<kol; i++) fscanf(FileIn,"%f,",&(sp->Xi[i]));


// читаємо масив Ak


for (i=0; i<kol; i++)


for (j=0; j<kol; j++) fscanf(FileIn,"%f,",&(sp->Ak[i][j]));


// читаємо n - кількість відрізків розбиття відрізку інтегрування


if(fscanf(FileIn,"%d,",&n)==NULL)


{ cout << strError; goto exit;};


down=sp->Xi[0];


up=sp->Xi[sp->kol-1]+(sp->Xi[sp->kol-1]-sp->Xi[sp->kol-2]);


I1=fSimps(down,up, n, sp);


I2=fSimps(down,up, 2*n, sp);


eps=(I2-I1)/15;


I=I2+eps;


// Вивід результатів


printf("I= %5.5fn",I);


printf("EPS= %5.5fn",eps);


fprintf(FileOut,"I= %5.5fn",I);


fprintf(FileOut,"EPS= %5.5fn",eps);


fclose(FileIn);


fclose(FileOut);


exit: cout << "n Press any key ...";


getch();}


Результат роботи програми ("data2_3ou.txt")


I= 1.32213


EPS= 0.00004


Завдання 3


Знайти розв’язок системи нелінійних рівнянь


,


Рішення.


Умову завдання перепишемо наступним чином


.


Приймаючи що і то коротко систему рівнянь можна записати так


.


Якщо відомо деяке наближення кореня системи рівнянь, то поправку можна знайти рішаючи систему


.


Розкладемо ліву частину рівняння по степеням малого вектору , обмежуючись лінійними членами


.


== – матриця похідних (матриця Якобі) ().


Складемо матрицю похідних (матрицю Якобі):



Якщо , то ,


де – матриця обернена до матриці Якобі.


Таким чином послідовне наближення кореня можна обчислити за формулою


або


.


Умовою закінчення ітераційного процесу наближення корення вибираємо умову


,


– евклідова відстань між двома послідовними наближеннями ;– число, що задає мінімальне наближення.


Для рішення систем нелінійних рівнянь за даним алгоритмом призначена програма


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.4fn",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.4fn f(x,y)= %4.3fn ",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= %lfn",data[0][0]);


for (int i=1;i<=CountF;i++)


printf("a%d = %lf , b%d = %lfn",i,data[i][0],i,data[i][1]);


cout << " Press any key ...";


getch();}


Результат роботи програми Work5.cpp




Сохранить в соц. сетях:
Обсуждение:
comments powered by Disqus

Название реферата: Числові методи

Слов:3469
Символов:40323
Размер:78.76 Кб.