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])
/* Функція обчислює коефіцієнти кубічного сплайну