Главная страница

Сплайны. Report_Рыбин_Хитева_cube_spline (1). Отчет по дисциплине Численные методы


Скачать 2.24 Mb.
НазваниеОтчет по дисциплине Численные методы
АнкорСплайны
Дата25.04.2022
Размер2.24 Mb.
Формат файлаdocx
Имя файлаReport_Рыбин_Хитева_cube_spline (1).docx
ТипОтчет
#494940

НИЖЕГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

ИНСТИТУТ РАДИОЭЛЕКТРОНИКИ И ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ

Кафедра «ПРИКЛАДНАЯ МАТЕМАТИКА»

ЛАБОРАТОРНАЯ РАБОТА №1

отчет

по дисциплине «Численные методы»

Одномерные и двумерные сплайны

Студент
__________ Рыбин А.В.

(Подпись) (Фамилия, И., О.)

_10-ПМ _________

(Группа) (Дата сдачи)
Студент
__________ Хитева Д.В.

(Подпись) (Фамилия, И., О.)

_10-ПМ _________

(Группа) (Дата сдачи)
Проверила
__________ Катаева Л.Ю.

(Подпись) (Фамилия, И.,О.)

Отчет защищён «____»_________2013г.
с оценкой____________________

Н. Новгород, 2013

Оглавление


Постановка задачи 3

Теоретический материал 4

Таблица идентификаторов 6

Реализация методов на С++ 6

Excel 21

Результат работы программ 23

Сводная таблица результатов 28

Вывод 29









Постановка задачи 3

Теоретический материал 4

Таблица идентификаторов 6

Реализация методов на С++ 6

Excel 21

Результат работы программ 23

Сводная таблица результатов 28

Вывод 29


Постановка задачи


Одним из наиболее широко используемых представлений кривых в компьютерном видении является В-сплайн представление. Важно различать сплайны и В-сплайны. В-сплайны являются полиномиальными функциями. Сплайны являются линейной комбинацией В-сплайнов. В литературе сплайны обычно определяются как различные виды степенной функции. Для вычислений более удобно определять сплайны рекурсивными функциями.

Некоторая функция f(x) задана на отрезке , разбитом на части , . Кубическим сплайном дефекта 1 называется функция , которая:

  • на каждом отрезке является многочленом степени не выше третьей;

  • имеет непрерывные первую и вторую производные на всём отрезке ;

  • в точках выполняется равенство , т. е. сплайн интерполирует функцию f в точках .

Для однозначного задания сплайна перечисленных условий недостаточно, для построения сплайна необходимо наложить какие-то дополнительные требования.

Естественным кубическим сплайном называется кубический сплайн, удовлетворяющий также граничным условиям вида:



Таким образом перед нами ставиться задача реализовать построение одномерного кубического сплайна. Так же мы ставим себе задаче реализовать построение двумерного кубического сплайна.

Теоретический материал


Обозначим: (1)

На каждом отрезке функция есть полином третьей степени , коэффициенты которого надо определить. Запишем для удобства в виде:



тогда



Условия непрерывности всех производных до второго порядка включительно записываются в виде






а условия интерполяции в виде



Отсюда получаем формулы для вычисления коэффициентов сплайна:








Формулы используемые при реализации одномерного и двумерного сплайна:

. (1)

. (2)
. (3)

при начальных условиях:

(4)

где N=n-1; n - количество точек.
Из заданных граничных условий, находим:
(5)
(6)



(7)

Полученные коэффициенты, это и есть нужные нам коэффициенты полиномов третьей степени.































Таблица идентификаторов


Таблица №1-Одномерный кубический сплайн

Тип переменной

Имя переменной

Хранимое значение

int

n

Кол-во точек

int

N

n-1

double

p[][]

Точки, координаты x и y

double

a,b

Начальные условия

double

m[]

коэффициенты

double

S[][]

Коэф-ты полученных уравнений для построения


Таблица №1-Двумерный кубический сплайн

Тип переменной

Имя переменной

Хранимое значение

int

kk

Кол-во столбцов

int

nn

Кол-во строк

double

p[][]

Точки, координаты

double

sy0,syn,sx0,sxn

Граничные условия

double

y0,yn,x0,xn

Начальные и конечные x и y

double

z[][]

Исходная матрица высот

double

zxy[][]

Итоговая расширенная матрица

Реализация методов на С++


Одномерный кубический сплайн

#include

#include
using namespace std;
int main()

{

double a,b; // init conditions

int n = 0; // кол-во точек

int N; //кол-во пар

int k = 0;
ifstream file;

file.open("in.txt");
while (file)

{

char ch = file.get();

if (ch == '\n')

n++;

}

file.close ();
double** p= new double*[n];

for(int i=0;i

{

p[i]=new double[2];

}
file.open("in.txt");

for (int i=0; i < n; ++i)

{

file>>p[i][0];

file>>p[i][1];

}

file.close ();
k=n;

cout<

cout<<"Enter initial conditions, S'("<


cin>>a>>b;

cout<<"S'("<


cout<<"points:"<

for(int i=0;i

{

cout<<"("<


}

/////////////////////

double* h = new double [n];

double* d = new double [n];

double* u = new double [n];

double* m = new double [n];

/////////////////////

for(int i=1;i

{

h[i-1]=p[i][0]-p[i-1][0]; // формула (1)

d[i-1]=(p[i][1]-p[i-1][1])/h[i-1]; // формула (2)

}

for(int i=1;i

{

u[i]=6*(d[i]-d[i-1]); // формула (3)

// cout<

}
N=n-1; //кол-во пар

double** M= new double*[n];

for(int i=1;i

{

M[i]=new double[n];

}

double** S= new double*[n];

for(int i=0;i

{

S[i]=new double[n];

}
// Формулы (4)

if(n==4)

{

for(int i=1;i

{

if(i==1)

{

M[1][1]=(3.0/2) * h[0] + 2 * h[1];

M[1][2]=h[1];

M[1][3]=u[1] - 3*(d[0]-a);

}

if(i>=2 && (N-2)==1)

{

M[i][1]=h[N-2];

M[i][2]=(2*h[N-2]+(3.0/2)*h[N-2]);

M[i][3]=u[N-1]-3*(b-d[N-1]);

}

}

}

if(n>4)

{

for(int i=1;i

{

if(i==1)

{

M[1][1]=(3.0/2) * h[0] + 2 * h[1];

M[1][2]=h[1];

for(int j=3;j

{

M[1][j]=0;

}

M[1][N]=u[1] - 3*(d[0]-a);

}

if((i>=2 && i<=N-1) && (N-2)!=1)

{

for(int j=1;j

{

if((i-1)!=j && i!=j && (i+1)!=j)

{

M[i][j]=0;

}

}

M[i][i-1]=h[i-1];

M[i][i]=2*(h[i-1]+h[i]);

M[i][i+1]=h[i];

M[i][N]=u[i];

}

if(i==N-1)

{

for(int j=1;j

{

if((N-2)!=j && (N-1)!=j)

{

M[i][j]=0;

}

}

M[i][N-2]=h[N-2];

M[i][N-1]=(2*h[N-2]+(3.0/2)*h[N-2]);

M[i][N]=u[N-1]-3*(b-d[N-1]);

}

}

}
for(int i=1;i

{

for(int j=1;j

{

cout<

}

cout<

}
n=n-2;

double temp;

//прямой ход

for(int i=1;i<=n;i++)

{

for(int j=i+1;j<=n;j++)

{

temp = (M[j][i])/(M[i][i]);

for(int k=1; k<=n+1;k++)

{

M[j][k]=M[j][k]-M[i][k]*temp;

}

}

}

//обратный ход

m[n] = M[n][n+1]/M[n][n];

//cout<
for(int i=n-1;i>0;i--)

{

temp=0;

for(int j=i+1;j<=n;j++)

{

temp = temp + M[i][j]*m[j];

}

temp = M[i][n+1] - temp;

m[i] = temp/M[i][i];

}

// Формулы (5) и (6)

m[0]=(3.0/h[0]) * (d[0] - a) - m[1]/2;

m[N]=(3.0/h[N-1]) * (b - d[N-1] ) - m[N-1]/2;
//ответ, вывод значения корней

for(int i=0;i

{

cout<<"m["<

}

// Формулы (7)

for(int i=1;i

{

S[i-1][0]=p[i-1][1];

S[i-1][1]=d[i-1]-(( h[i-1] * (2*m[i-1]+m[i]))/6);

S[i-1][2]=m[i-1]/2;

S[i-1][3]=(m[i]-m[i-1])/(6*h[i-1]);

cout<<"S["<

}
fstream outfile;

outfile.open("spline.txt", ios::out); //ios::app ios::in

outfile<

for(int i=0;i

{

outfile<<(p[i][0])<<" ";

}

outfile<<"\n";

for(int i=1;i

{

outfile<

}

outfile.close();
system("pause");

return 0;

}



Двумерный кубический сплайн

#include

#include

#include
using namespace std;
int main()

{

int kk = 0; // кол-во столбцов

int nn = 0; // количество строк

double y0,yn,x0,xn; // начальные x и y, и конечные x и y.

double sx0,sxn; // начальные условия для x

double sy0,syn; // начальные условия для y

int N;

int n = 0;

int jj;

int mm;

int iy;

int jy;

double hx;

double hy;
ifstream file;

file.open("in1.txt");
while (file)

{

char ch = file.get();

if(ch=='\t')

{

kk++;

}

if (ch == '\n' )

{

nn++;

}

}

file.close ();
double** z= new double*[nn];

for(int i=0;i

{

z[i]=new double[kk];

}
file.open("in1.txt");

for (int i=0; i < nn; ++i)

{

for(int j=0;j

{

file>>z[i][j];

}

}

file.close ();
for (int i=0; i < nn; ++i)

{

for(int j=0;j

{

cout<

}

cout<

}
double** p= new double*[nn];

for(int i=0;i

{

p[i]=new double[2];

}
double** zx= new double*[5*nn];

for(int i=0;i<5*nn;++i)

{

zx[i]=new double[kk];

}

double* x = new double [nn];

double* x_point = new double [5*nn];

double* h = new double [nn];

double* d = new double [nn];

double* u = new double [nn];

double* m = new double [nn];

double** M= new double*[nn];

for(int i=1;i

{

M[i]=new double[nn];

}
double** S= new double*[nn];

for(int i=0;i

{

S[i]=new double[4];

}
//cout<

//cout<

// x0=1;

//xn=8;

// y0=1;

// yn=10;

// sx0=0.2;

// sxn=-1;

// sy0=0.3;

// syn=-2;
cout<<"Enter x0 and xn, y0 and yn"<

cin>>x0>>xn;

cin>>y0>>yn;
cout<<"Enter initial conditions for x , S'("<

cin>>sx0>>sxn;
cout<<"Enter initial conditions for y , S'("<

cin>>sy0>>syn;

///////////////// x
hx=(xn-x0)/(nn-1);

//cout<
x[0]=x0;

double temp=x0;

for(int i=1;i

{

temp=temp+hx;

x[i]=temp;

}
for(int imt=0;imt

{
for(int i=0;i

{

p[i][0]=x[i];

p[i][1]=z[i][imt];

// cout<<"("<


}
for(int i=1;i

{

h[i-1]=p[i][0]-p[i-1][0]; // формула (1)

d[i-1]=(p[i][1]-p[i-1][1])/h[i-1]; // формула (2)

//cout<

//cout<

}
for(int i=1;i

{

u[i]=6*(d[i]-d[i-1]); // формула (3)

//cout<

}

//

N=nn-1;

// формула (4)

if(nn==4)

{

for(int i=1;i

{

if(i==1)

{

M[1][1]=(3.0/2) * h[0] + 2 * h[1];

M[1][2]=h[1];

M[1][3]=u[1] - 3*(d[0]-sx0);

}

if(i>=2 && (N-2)==1)

{

M[i][1]=h[N-2];

M[i][2]=(2*h[N-2]+(3.0/2)*h[N-2]);

M[i][3]=u[N-1]-3*(sxn-d[N-1]);

}

}

}
if(nn>4)

{

for(int i=1;i

{

if(i==1)

{

M[1][1]=(3.0/2) * h[0] + 2 * h[1];

M[1][2]=h[1];

for(int j=3;j

{

M[1][j]=0;

}

M[1][N]=u[1] - 3*(d[0]-sx0);

}

if((i>=2 && i<=N-1) && (N-2)!=1)

{

for(int j=1;j

{

if((i-1)!=j && i!=j && (i+1)!=j)

{

M[i][j]=0;

}

}

M[i][i-1]=h[i-1];

M[i][i]=2*(h[i-1]+h[i]);

M[i][i+1]=h[i];

M[i][N]=u[i];

}

if(i==N-1)

{

for(int j=1;j

{

if((N-2)!=j && (N-1)!=j)

{

M[i][j]=0;

}

}

M[i][N-2]=h[N-2];

M[i][N-1]=(2*h[N-2]+(3.0/2)*h[N-2]);

M[i][N]=u[N-1]-3*(sxn-d[N-1]);

}

}

}

n=nn-2;

//double temp;

temp=0;

//прямой ход

for(int i=1;i<=n;i++)

{

for(int j=i+1;j<=n;j++)

{

temp = (M[j][i])/(M[i][i]);

for(int k=1; k<=n+1;k++)

{

M[j][k]=M[j][k]-M[i][k]*temp;

}

}

}

//обратный ход

m[n] = M[n][n+1]/M[n][n];

//cout<
for(int i=n-1;i>0;i--)

{

temp=0;

for(int j=i+1;j<=n;j++)

{

temp = temp + M[i][j]*m[j];

}

temp = M[i][n+1] - temp;

m[i] = temp/M[i][i];

}

// формула (5) и (6)

m[0]=(3.0/h[0]) * (d[0] - sx0) - m[1]/2;

m[N]=(3.0/h[N-1]) * (sxn - d[N-1]) - m[N-1]/2;

for(int i=1;iФормула (7)

{

S[i-1][0]=p[i-1][1];

S[i-1][1]=d[i-1]-((h[i-1]*(2*m[i-1]+m[i]))/6);

S[i-1][2]=m[i-1]/2;

S[i-1][3]=(m[i]-m[i-1])/(6*h[i-1]);

}

jj=(xn-x0)/(hx/5);

x_point[0]=x0;

temp=x0;

for(int i=1;i<=jj;++i)

{

temp=temp+(hx/5);

x_point[i]=temp;

//cout<

}
mm=0;

for(int i=1;i

{

for(int j=0;j<5;j++)

{

zx[mm][imt]=( S[i-1][3] )*pow((x_point[mm]-p[i-1][0]),3) +(S[i-1][2])*pow((x_point[mm]-p[i-1][0]),2) +( S[i-1][1] )*( x_point[mm] - p[i-1][0]) + ( S[i-1][0] );

// cout<

mm=mm+1;

}

}

zx[mm][imt]=z[nn-1][imt];

}
///////////////////

double** py= new double*[kk];

for(int i=0;i

{

py[i]=new double[2];

}
double* y= new double [kk];

double* y_point = new double [5*kk];
double* h1 = new double [kk];

double* d1 = new double [kk];

double* u1 = new double [kk];

double* m1 = new double [kk];

double** M1= new double*[kk];

for(int i=1;i

{

M1[i]=new double[kk];

}
double** S1= new double*[kk];

for(int i=0;i

{

S1[i]=new double[4];

}
double** zxy= new double*[5*nn];

for(int i=0;i<5*nn;++i)

{

zxy[i]=new double[5*kk];

}

//////////// по y
hy=(yn-y0)/(kk-1);

//cout<
y[0]=y0;

temp=y0;
for(int i=1;i

{

temp=temp+hy;

y[i]=temp;

//cout<

}
for(int jmt=0;jmt<=jj;jmt++)

{
for(int i=0;i

{

py[i][0]=y[i];

py[i][1]=zx[jmt][i]; //jmt

// cout<<"("<


}

for(int i=1;i

{

h1[i-1]=py[i][0]-py[i-1][0];

d1[i-1]=(py[i][1]-py[i-1][1])/h1[i-1];

//cout<

//cout<

}

for(int i=1;i

{

u1[i]=6*(d1[i]-d1[i-1]);

//cout<

}
N=kk-1;

//

if(kk==4)

{

for(int i=1;i

{

if(i==1)

{

M1[1][1]=(3.0/2) * h1[0] + 2 * h1[1];

M1[1][2]=h1[1];

M1[1][3]=u1[1] - 3*(d1[0]-sy0);

}

if(i>=2 && (N-2)==1)

{

M1[i][1]=h1[N-2];

M1[i][2]=(2*h1[N-2]+(3.0/2)*h1[N-2]);

M1[i][3]=u1[N-1]-3*(syn-d1[N-1]);

}

}

}
if(kk>4)

{

for(int i=1;i

{

if(i==1)

{

M1[1][1]=(3.0/2) * h1[0] + 2 * h1[1];

M1[1][2]=h1[1];

for(int j=3;j

{

M1[1][j]=0;

}

M1[1][N]=u1[1] - 3*(d1[0]-sy0);

}

if((i>=2 && i<=N-1) && (N-2)!=1)

{

for(int j=1;j

{

if((i-1)!=j && i!=j && (i+1)!=j)

{

M1[i][j]=0;

}

}

M1[i][i-1]=h1[i-1];

M1[i][i]=2*(h1[i-1]+h1[i]);

M1[i][i+1]=h1[i];

M1[i][N]=u1[i];

}

if(i==N-1)

{

for(int j=1;j

{

if((N-2)!=j && (N-1)!=j)

{

M1[i][j]=0;

}

}

M1[i][N-2]=h1[N-2];

M1[i][N-1]=(2*h1[N-2]+(3.0/2)*h1[N-2]);

M1[i][N]=u1[N-1]-3*(syn-d1[N-1]);

}

}

}
n=kk-2;

//double temp;

temp=0;

//прямой ход

for(int i=1;i<=n;i++)

{

for(int j=i+1;j<=n;j++)

{

temp = (M1[j][i])/(M1[i][i]);

for(int k=1; k<=n+1;k++)

{

M1[j][k]=M1[j][k]-M1[i][k]*temp;

}

}

}

//обратный ход

m1[n] = M1[n][n+1]/M1[n][n];

//cout<
for(int i=n-1;i>0;i--)

{

temp=0;

for(int j=i+1;j<=n;j++)

{

temp = temp + M1[i][j]*m1[j];

}

temp = M1[i][n+1] - temp;

m1[i] = temp/M1[i][i];

}
m1[0]=(3.0/h1[0]) * (d1[0] - sy0) - m1[1]/2;

m1[N]=(3.0/h1[N-1]) * (syn - d1[N-1]) - m1[N-1]/2;
for(int i=1;i

{

S1[i-1][0]=py[i-1][1];

S1[i-1][1]=d1[i-1]-((h1[i-1]*(2*m1[i-1]+m1[i]))/6);

S1[i-1][2]=m1[i-1]/2;

S1[i-1][3]=(m1[i]-m1[i-1])/(6*h1[i-1]);

}


jy=(yn-y0)/(hy/5);

y_point[0]=y0;

temp=y0;

for(int i=1;i<=jy;++i)

{

temp=temp+(hy/5);

y_point[i]=temp;

}
iy=0;

for(int i=1;i

{

for(int j=0;j<5;j++)

{

zxy[jmt][iy]=( S1[i-1][3] )*pow((y_point[iy]-py[i-1][0]),3) +(S1[i-1][2])*pow((y_point[iy]-py[i-1][0]),2) +( S1[i-1][1] )*( y_point[iy] - py[i-1][0]) + ( S1[i-1][0] );

iy=iy+1;

}

}

zxy[jmt][iy]=zx[jmt][kk-1]; //jmt

}
fstream outfile;

outfile.open("spline_2D.txt", ios::out | ios::in); //ios::app ios::in ios::app |

outfile<< jj <<" "<< jy <<" \n";

for(int i=0;i<=jj;i++)

{

outfile<

}

outfile<<"\n";

for(int i=0;i<=jy;i++)

{

outfile<

}

outfile<<"\n";

for(int i=0;i<=jj;i++)

{

for(int j=0;j<=jy;j++)

{

outfile<

}

outfile<<"\n";

}

outfile<< nn <<" "<< kk <<" \n";

for(int i=0;i

{

outfile<

}

outfile<<"\n";

for(int i=0;i

{

outfile<

}

outfile<<"\n";
for(int i=0;i

{

for(int j=0;j

{

outfile<

}

outfile<<"\n";

}
outfile.close();

cout<<"ok"<
/////////////

delete []d;

delete []m;

delete []u;

delete []h;

delete []x;
for(int i=0;i

delete [] z[i];

delete [] z;
for(int i=0;i

delete [] S[i];

delete [] S;
for(int i=0;i

delete [] zx[i];

delete [] zx;

///////////////

system("pause");

return 0;

}


Excel
Реализация нахождения коэффициентов одномерного кубического сплайна средствами Excel (Хитева Д.В.).

На рис.1 представлена реализация нахождения коэффициентов кубического сплайна для конкретного примера по четырем точкам. В строках 2-3, а так же в клетках H5, H6, I5, I6 представлены исходные данные задачи. В строках 5-13 производятся вычисления. В строках 15-18 записан ответ.



Рис.1 Кубический сплайн
На рис.1.1 представлена реализация указанного выше задания в формулах.



Рис.1.1 Кубический сплайн в формулах


Реализация нахождения коэффициентов одномерного кубического сплайна средствами Excel (Рыбин А.В.).

На рис.2 представлена реализация нахождения коэффициентов кубического сплайна для конкретного примера по четырем точкам. В строках 2-3, а так же в клетках H5, H6, I5, I6 представлены исходные данные задачи. В строках 5-13 производятся вычисления. В строках 15-18 записан ответ.



Рис.2 Кубический сплайн
На рис.2.1 представлена реализация указанного выше задания в формулах.



Рис.2.1 Кубический сплайн в формулах

Результат работы программ


Реализация на С++ (Хитева Д.В.).

На рис.3 представлена иллюстрация к реализации одномерного кубического сплайна.



Рис.-3 Одномерный кубический сплайн
На рис.4 представлена иллюстрация к реализации двумерного кубического сплайна.



Рис.-4 Двумерный кубический сплайн

Реализация в пакете SciLab (Хитева Д.В.).

На рис.5 представлена иллюстрация построения одномерного кубического сплайна.



Рис.-5 Одномерный кубический сплайн
На рис.6 представлена иллюстрация построения двумерного кубического сплайна.



Рис.-6 Двумерного кубический сплайн
Реализация на С++ (Рыбин А.В).

На рис.7 представлена иллюстрация к реализации одномерного кубического сплайна.



Рис.-7 Одномерный кубический сплайн
На рис.8 представлена иллюстрация к реализации двумерного кубического сплайна.



Рис.-8 Одномерный кубический сплайн


Реализация в пакете SciLab (Рыбин А.В).

На рис.9 представлена иллюстрация построения одномерного кубического сплайна.



Рис.-9 Одномерный кубический сплайн

На рис.10 представлена иллюстрация построения двумерного кубического сплайна.



Рис.-10 Двумерный кубический сплайн


Сводная таблица результатов



Таблица №1. Сводная таблица результатов для одномерного сплайна (Хитева Д.В.)




S00

S01

S02

S03

S10

S11

S12

S13

S20

S21

S22

S23

Excel

1

1

0,4

0,6

3

3,6

2,2

-2,8

6

-0,4

-6,2

2,6

C++

1

1

0,4

0,6

3

3,6

2,2

-2,8

6

-0,4

-6,2

2,6


Таблица №2. Сводная таблица результатов для одномерного сплайна (Рыбин А.В.)




S00

S01

S02

S03

S10

S11

S12

S13

S20

S21

S22

S23

Excel

1

3

-2,73

2,73

4

5,73

5,47

-6,2

9

-1,93

-13,1

8,07

C++

1

3

-2,73

2,73

4

5,73

5,47

-6,2

9

-1,93

-13,1

8,07


Вывод


Форма сплайна задается кубическим полиномом.

В общем случае математический сплайн - это кусочный полином степени k с непрерывной производной степени k-1 в точках соединения сегментов. Так, например, кубический сплайн имеет в точках соединения непрерывность второго порядка. Кусочные сплайны из многочленов невысокого порядка очень удобны для интерполяции кривых, так как они не требуют больших вычислительных затрат и не вызывают численных отклонений, свойственных многочленам высокого порядка. По аналогии с физическими сплайнами обычно используется серия кубических сегментов, причем каждый сегмент проходит через две точки. Кубический сплайн удобен еще и тем, что это кривая наименьшего порядка, допускающая точки перегиба и изгиб в пространстве.


написать администратору сайта