Числові методи (стр. 3 из 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,"%d\n",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.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 << "\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 - вказівник на функцію, яка обчислює по