Числові методи (стр. 2 из 4)
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])
/* Функція обчислює коефіцієнти кубічного сплайну