НИЖЕГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ
ИНСТИТУТ РАДИОЭЛЕКТРОНИКИ И ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ
Кафедра «ПРИКЛАДНАЯ МАТЕМАТИКА»
ЛАБОРАТОРНАЯ РАБОТА №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 в точках соединения сегментов. Так, например, кубический сплайн имеет в точках соединения непрерывность второго порядка. Кусочные сплайны из многочленов невысокого порядка очень удобны для интерполяции кривых, так как они не требуют больших вычислительных затрат и не вызывают численных отклонений, свойственных многочленам высокого порядка. По аналогии с физическими сплайнами обычно используется серия кубических сегментов, причем каждый сегмент проходит через две точки. Кубический сплайн удобен еще и тем, что это кривая наименьшего порядка, допускающая точки перегиба и изгиб в пространстве. |