МІНІСТЕРСТВООСВІТИ УКРАЇНИ
ЧЕРНІВЕЦЬКИЙДЕРЖАВНИЙ УНІВЕРСИТЕТ
ІМ. Ю.ФЕДЬКОВИЧА
КОНТРОЛЬНАРОБОТА
здисципліни " Числові методи "
Варіант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
#include
#include
const int nMax=5;// максимальна кількість рівнянь
const floatZERO=.0000001;
int fGaus(floatA[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
{aMax=A[k][k]; l=k;
for (i=k+1;i
if(fabs(A[i][k])>fabs(aMax))
{aMax=A[i][k];
l=i;}
// якщо модульголовного елементу aMax менший за програмний 0 (ZERO)
if ( fabs(aMax)
// якщо потрібно,міняємо місцями рівняння Pk i Pl
if ( l!=k)
{for( j=0;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
B[k]/=aMax;
// обчислюємо коефіцієнтиA[i][j] та вільні члени решти рівнянь
for (i=k+1;i
{t=A[i][k];B[i]-=t*B[k];
for (j=0; j
} // for (k)
// Зворотній хід
for ( k=n-1;k>=0; k--)
{X[k]=0;
for (l=k+1;l
X[k]=B[k]-X[k];}
return 0;
} // fGaus()
void main()
{floatA[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
goto exit;}
FileOut=fopen(«data_out.txt»,«w»);// відкриваємо файл для запису
if(FileOut==NULL)
{cout
goto exit;}
if(fscanf(FileIn,"%d",&n)==NULL)
{ cout
for (i=0; i
for(j=0; j
fscanf(FileIn,"%f",&(A[i][j]));
for (i=0;i
if(fscanf(FileIn,"%f",&(B[i]))==NULL)
{ cout
if(fGaus(A,B,n,X)!=0)
{ cout
// Вивідрезультатів
for (i=0; i
{printf("x[%d]= %f ",i+1,X[i]);
fprintf(FileOut,"x[%d]= %f ",i+1,X[i]);}
fclose(FileIn);
fclose(FileOut);
exit: cout
getch();}
Результат роботипрограми:
x[1]= 3.017808 x[2]=0.356946 x[3]= -0.302131
б) Знайдемообернену матрицю />.
0-й крок.
А /> Е />
1-й крок.
/>
/> />
/>
/> />
/> ; />
/> />
2-й крок.
/>
/> />
/> ; />
/> />
3-й крок.
/>; /> ;/>
/> />
/>
/>
/>
/>.
Даний алгоритмрішення системи лінійних рівнянь реалізований в програмі Work1_2.
//------------------------------------------------------------
// Work1_2.cpp
//------------------------------------------------------------
// «Числовіметоди»
// Завдання 1
// Рішеннясистеми лінійних рівнянь методом Гауса-Жордана
#include
#include
#include
const int nMax=5;// максимальна кількість рівнянь
const floatZERO=.0000001;
intfGausJordan(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
for (j=0; j
Ainv[i][j] =(i==j)? 1.: 0.;
for (k=0; k
{// знаходимо махпо модулю елемент
aMax=A[k][k];l=k;
for (i=k+1;i
if(fabs(A[i][k])>fabs(aMax))
{ aMax=A[i][k]; l=i;}
// якщо модульголовного елементу aMax менший за програмний 0 (ZERO)
if (fabs(aMax)
// якщо потрібно,міняємо місцями рівняння Pk i Pl
if ( l!=k)
for( j=0; 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
// обчислюємоелементи решти рядків
for (i=0; i
if( i!=k )
{t=A[i][k];
for (j=0; j
{A[i][j]-=t*A[k][j];
Ainv[i][j]-=t*Ainv[k][j];}}}
return 0;
} //fGausJordana()
void fDobMatr(intn, float A[nMax][nMax], float B[nMax],float X[nMax])
// функціязнаходить добуток матриці А на вектор В і результат повертає в
// векторі Х
{int i,j;
float summa;
for (i=0; i
{summa=0;
for (j=0; j
{summa+=A[i][j]*B[j];
X[i]=summa;}}
} // fDobMatr
void main()
{floatA[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
goto exit;}
FileOut=fopen(«data_out.txt»,«w»);// відкриваємо файл для запису
if(FileOut==NULL)
{cout
goto exit;}
if(fscanf(FileIn,"%d",&n)==NULL)
{ cout
for (i=0; i
for(j=0; j
fscanf(FileIn,"%f",&(A[i][j]));
for (i=0;i
if(fscanf(FileIn,"%f",&(B[i]))==NULL)
{ cout
if(fGausJordan(n,A,Ainv)!=0)
{ cout
fDobMatr(n,Ainv,B,X);
// Вивідрезультатів
for (i=0; i
{printf("x[%d]= %f ",i+1,X[i]);
fprintf(FileOut,"x[%d]= %f ",i+1,X[i]);}
fclose(FileIn);
fclose(FileOut);
exit: cout
getch();}
Результат роботипрограми:
x[1]= 3.017808 x[2]=0.356946 x[3]= -0.302131
Завдання 2
Задана задачаКоші
/>, />
а) Знайтирозв’язок /> в табличній формі методомРунге-Кутта:
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>
/>, />,/> />.
б) Інтерполюватицю функцію кубічним сплайном. Систему рівнянь для моментів кубічного сплайнурозв’язати методом прогонки. Вибрати крайові умови для кубічного сплайну у вигляді
/>.
в) Використовуючикубічний сплайн з пункту б) обчислити /> методомСімпсона .
Взяти /> (/>– кількість відрізківрозбиття).
Рішення.
а) МетодРунге-Кутта
Розрахунок будемопроводити за наступними формулами :
/>;
/>;
/>;
/>;
/>;
/>.
Цей алгоритмреалізовується в програмі Work2_1.
//------------------------------------------------------------
// Work2_1.cpp
//------------------------------------------------------------
// «Числовіметоди»
// Завдання 2
// Рішення задачіКоші методом Рунге-Кутта
#include
#include
#include
typedef float(*pfunc)(float,float); // pfunc — вказівник на функцію
const int nMax=5;// максимальна кількість відрізків розбиття
voidfRunge_Kutta(pfunc f, float x0, float y0,float h, int n, float Y[nMax])
/* Функціязнаходить табличне значення функції методом Рунге-Кутта
Вхідні дані:
f — функціяf(x,y)
x0,y0 — початковаточка;
h — крок;
n- кількістьточок розбиття;
Вихідні дані:
Y- вектор значеньфункції*/
{floatk1,k2,k3,k4,x; // максимальний елемент, тимчасова змінна
int i;
x=x0; Y[0]=y0;
for (i=0;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;}}
floatMyfunc(float x,float y)
{returnlog10(cos(x+y)*cos(x+y)+2)/log10(5);}
void main()
{floatY[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
goto exit;}
FileOut=fopen(«data2_out.txt»,«w»);// відкриваємо файл для запису
if(FileOut==NULL)
{cout
goto exit;}
FileOut2=fopen(«data2_2in.txt»,«w»);// відкриваємо файл для запису
if(FileOut==NULL)
{cout
goto exit;}
if(fscanf(FileIn,"%d%f%f%f,",&n,&h,&x0,&y0)==NULL)
{ cout
fRunge_Kutta(Myfunc,x0,y0,h,n,Y);
// Вивідрезультатів
for (i=0; 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
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
#include
#include
const int nMax=4;// максимальна кількість відрізків розбиття
const floatx0=0.;// початкова точка сітки
const floath=0.1;// крок розбиття
// векториматриці А
float a[]={0.,0.5, 0.5};
float b[]={2.,2., 2.};
float c[]={0.5,0.5, 0.};
//voidfMetodProgonku( int n,float a[nMax],float b[nMax],float c[nMax],float d[nMax],float M[nMax+1])
/* Функція знаходитьмоменти кубічного сплайну методом прогонки
Вхідні дані:
a,b,c -векториматриці А ;
d — векторвільних членів;
n- степіньматриці А;
Вихідні дані:
М- вектормоментів кубічного сплайну.*/
{floatk[nMax],fi[nMax];
int i;
// прямий хід
for (i=0; 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( intn,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
{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()
{floatY[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
goto exit; }
FileOut=fopen(«data2_2ou.txt»,«w»);// відкриваємо файл для запису
if(FileOut==NULL)
{ cout
goto exit; }
FileOut2=fopen(«data2_3in.txt»,«w»);// відкриваємо файл для запису
if(FileOut2==NULL)
{ cout
goto exit; }
// читаємо векторY
for (i=0;i
if(fscanf(FileIn,"%f,",&(Y[i]))==NULL)
{ cout
// обчислюємовектор d
for (i=1; i
//fMetodProgonku(n-1,a,b,c,d,M);
fSplain(n,h,Y,M,Ak);
// Вивідрезультатів в тому числі і для наступного завдання
fprintf(FileOut2,"%d\n",n);// n — кількість відрізків
// координатиточок сітки по Х
for(floatxi=x0,i=0; i
fprintf(FileOut2,"\n");
for (i=0; i
{for (j=0;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
fprintf(FileOut,"\n");
fprintf(FileOut2,"\n");}
fclose(FileIn);
fclose(FileOut);
exit: cout
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
#include
#include
#include
// визначеннясплайнового класу
class Tsplain
{public:
int kol; //кількість рівнянь (відрізків розбиття)
float ** Ak; //масив коефіцієнтів
float * Xi; //вектор початків відрізків
float vol(floatx); // функція повертає значення сплайну в точці х
Tsplain(int k);// constructor};
Tsplain::Tsplain(intk)
{kol=k;
Xi=newfloat[kol];
Ak=newfloat*[kol];
for(int i=0;i
floatTsplain::vol(float x)
{float s=0.;
int i,t;
// шукаємовідрізок t де знаходиться точка х
for (i=0;i=Xi[i]) { t=i; break; }
s=Ak[t][0];
for (i=1;i
s+=Ak[t][i]*pow((x-Xi[t]),i);
return s;}
floatfSimps(float down,float up, int n, Tsplain *spl)
/* Функція обчислюєінтеграл методом Сімпсона з використанням кубічного сплайну
Вхідні дані:
down,up -границіінтегрування ;
n- числовідрізків, на яке розбиваєтьться відрізок інтегрування ;
spl — вказівникна об’єкт класу Tsplain ( кубічний сплайн );
Вихідні дані:
функція повертаєзнайдене значення інтегралу.*/
{float s=0;
floath=(up-down)/(float)n;
int i;
s=spl->vol(down)+spl->vol(up-h);
for (i=2; i
s+=2*(spl->vol(down+i*h));
for (i=1; i
s+=4*(spl->vol(down+i*h));
return s*h;}
void main()
{int kol; //кількість рівняннь кубічного сплайну
float down,up;
floatI1,I2,I,eps;
int n,i,j;
char *strError="\nError of file !";
FILE*FileIn,*FileOut;
FileIn=fopen(«data2_3in.txt»,«r»);// відкриваємо файл для читання
if (FileIn==NULL)
{ cout
goto exit; }
FileOut=fopen(«data2_3ou.txt»,«w»);// відкриваємо файл для запису
if(FileOut==NULL)
{ cout
goto exit; }
// читаємо kol
if(fscanf(FileIn,"%d,",&kol)==NULL)
{ cout
Tsplain *sp;
sp=newTsplain(kol);
// читаємо векторXi
for(i=0;iXi[i]));
// читаємо масивAk
for (i=0;i
for (j=0;jAk[i][j]));
// читаємо n — кількість відрізків розбиття відрізку інтегрування
if(fscanf(FileIn,"%d,",&n)==NULL)
{ cout
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.5f\n»,I);
printf(«EPS=%5.5f\n»,eps);
fprintf(FileOut,«I=%5.5f\n»,I);
fprintf(FileOut,«EPS=%5.5f\n»,eps);
fclose(FileIn);
fclose(FileOut);
exit: cout
getch();}
Результат роботипрограми («data2_3ou.txt»)
I= 1.32213
EPS= 0.00004
Завдання 3
Знайти розв’язоксистеми нелінійних рівнянь
/> , />
Рішення.
Умову завданняперепишемо наступним чином
/> .
Приймаючи що /> і /> то коротко систему рівняньможна записати так
/>.
Якщо відомо деякенаближення /> кореня /> системи рівнянь, топоправку /> можна знайти рішаючисистему
/>.
Розкладемо лівучастину рівняння по степеням малого вектору />,обмежуючись лінійними членами
/> /> /> .
/>=/>=/> – матриця похідних(матриця Якобі) (/>).
Складемо матрицюпохідних (матрицю Якобі):
/>
Якщо /> , то /> ,
де /> – матриця обернена доматриці Якобі.
Таким чиномпослідовне наближення кореня можна обчислити за формулою
/> або
/>.
Умовою закінченняітераційного процесу наближення корення вибираємо умову
/>,
/>– евклідова відстань міждвома послідовними наближеннями ;/>– число,що задає мінімальне наближення.
Для рішеннясистем нелінійних рівнянь за даним алгоритмом призначена програма
Work3.cpp
//------------------------------------------------------------
// Work3.cpp
//------------------------------------------------------------
// «Числовіметоди»
// Завдання 3
// Розв’язуваннясистеми нелінійних рівнянь методом Ньютона
#include
#include
#include
#include
#include«matrix.h»
const int N=2; //степінь матриці Якобі (кількість рівнянь)
typedef void(*funcJ) (float[N], float[N][N]);
voidfJakobi(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(floatX[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];}
//intNelinSist(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;
floatW[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
printf(«X=%5.4f Y= %5.4f\n»,X[0],X[1]);
cout
getch();}
Результат роботипрограми:
X= 0.1477 Y=1.0214
Завдання 4
Знайти точку мінімумута мінімальне значення функції
/> , />
методом Ньютона.
Рішення.
/>;
Матриця Гессе
/>.
Ітераційнийпроцес послідовного наближення мінімуму функції буде таким:
/>,
де />– матриця обернена доматриці Гессе.
Для закінченняітераційного процесу використаємо умову
/> або
/>.
Для пошуку мінімумуфункції за методом Ньютона призначена програма Work4.cpp
//------------------------------------------------------------
// matrix.h
//-----------------------------------------------------------
const int nMax=2;// кількість рівнянь
const floatZERO=.00000001;
int invMatr(intn,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
for (j=0; j
Ainv[i][j] =(i==j)? 1.: 0.;
for (k=0; k
{// знаходимо махпо модулю елемент
aMax=A[k][k];l=k;
for (i=k+1;i
if(fabs(A[i][k])>fabs(aMax))
{ aMax=A[i][k]; l=i;}
// якщо модульголовного елементу aMax менший за програмний 0 (ZERO)
if (fabs(aMax)
// якщо потрібно,міняємо місцями рівняння Pk i Pl
if ( l!=k)
for( j=0; 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
// обчислюємоелементи решти рядків
for (i=0; i
if( i!=k )
{t=A[i][k];
for (j=0; j
{A[i][j]-=t*A[k][j];
Ainv[i][j]-=t*Ainv[k][j];}}}
return 0;}
void DobMatr(intn, float A[nMax][nMax], float B[nMax],float X[nMax])
// функціязнаходить добуток матриці А на вектор В і результат повертає в
// векторі Х
{int i,j;
float summa;
for (i=0; i
{summa=0;
for (j=0; 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; //степінь матриці Гессе
floatmyFunc(float x[N])
{ returnexp(-x[1])-pow(x[1]+x[0]*x[0],2); }
typedef void(*funcH) (float[N], float[N][N]);
void fHesse(floatX[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(floatX[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(floatX[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;
floatHesse[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
printf(«X=%5.5f Y= %5.4f\n f(x,y)= %4.3f\n »,X[0],X[1],myFunc(X));
cout
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
const doublePi=3.141592653;
// функціяповертає і-й вузол квадратурної формули, 2N+1-кілікість вузлів
inline doubleFuncXi(int N, int i) {return -Pi+(2*Pi*i)/(2*N+1);}
typedef double(*Func)(double); // опис типу вказівника на функцію
char Fourier(FuncF_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
Fourier=2 — якщоCountK
{doublea,b,sumA,sumB;
int i,k;
if (CountN
if (CountK
// обчислення а0
sumA=0;
for (i=0; i
a=1./(2*CountN+1)*sumA;
Arr[0][0]=a;
// обчисленнякоефіцієнтів аk,bk
for (k=1;k
{sumA=sumB=0;
for (i=0;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
#include
#include
double f(doublex)
// функціяповертає значення функції f(x)
{returnsqrt(Pi*Pi*x*x+1);}
const int N=20;// константа, яка визначає розбиття відрізка [-п; п]
// на рівнічастини
const intCountF=15; // кількість пар коефіцієнтів ряду
void main()
{double **data;
data = new double*[CountF+1];
for ( int i=0;i
if(Fourier(f,N,CountF,data) != 0)
{cout
return;}
// Вивідрезультатів
printf(«a0=%lf\n»,data[0][0]);
for (inti=1;i
printf(«a%d= %lf, b%d = %lf\n»,i,data[i][0],i,data[i][1]);
cout
getch();}
Результат роботипрограми Work5.cpp
/>