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,"%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

Знайти розв'язок системи нелінійних рівнянь

,

Страницы: 1, 2, 3