МІНІСТЕРСТВО ОСВІТИ УКРАЇНИ
ЧЕРНІВЕЦЬКИЙ ДЕРЖАВНИЙ УНІВЕРСИТЕТ
ІМ. Ю. ФЕДЬКОВИЧА
КОНТРОЛЬНА РОБОТА
з дисципліни " Числові методи "
Варіант 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(
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