НИЖЕГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ
УНИВЕРСИТЕТ
ИНСТИТУТ РАДИОЭЛЕКТРОНИКИ И ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ
Кафедра «ПРИКЛАДНАЯ МАТЕМАТИКА»
Лабораторная работа №1
отчет
по дисциплине «Численные методы»
Одномерные и двумерные сплайны
Студент
__________ Рыбин А.В.
(Подпись) (Фамилия, И., О.)
_10-ПМ _________
(Группа) (Дата сдачи)
Студент
__________ Хитева Д.В.
(Подпись) (Фамилия, И., О.)
_10-ПМ _________
(Группа) (Дата сдачи)
Проверила
__________ Катаева Л.Ю.
(Подпись) (Фамилия, И.,О.)
Отчет защищён «____»_________2013г.
с оценкой____________________
Н. Новгород, 2013
Оглавление
Постановка задачи 2
Теоретический материал 4
Таблица идентификаторов 6
Реализация методов на С++ 7
Excel 21
Результат работы программ 23
Сводная таблица результатов 28
Вывод 29
Постановка задачи 3
Теоретический материал 4
Таблица идентификаторов 6
Реализация методов на С++ 7
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 <iostream> #include <fstream>
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<n;++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<<n<<endl; cout<<"Enter initial conditions, S'("<<p[0][0]<<") and S'("<<p[n-1][0]<<")"<<endl; cin>>a>>b; cout<<"S'("<<p[0][0]<<")="<< a <<" and S'("<<p[n-1][0]<<")="<< b <<endl<<endl;
cout<<"points:"<<endl; for(int i=0;i<n;i++) { cout<<"("<<p[i][0]<<","<<p[i][1]<<")"<<endl; } ///////////////////// 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<n;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<n-1;i++) { u[i]=6*(d[i]-d[i-1]); // формула (3) // cout<<u[i]<<endl; }
N=n-1; //кол-во пар
double** M= new double*[n]; for(int i=1;i<N;++i) { M[i]=new double[n]; }
double** S= new double*[n]; for(int i=0;i<k;++i) { S[i]=new double[n]; }
// Формулы (4) if(n==4) { for(int i=1;i<N;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<N;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<n;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<n;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<n;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<N;i++) { for(int j=1;j<n;j++) { cout<<M[i][j]<<" "; } cout<<endl; }
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<<n<<endl;
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<k;i++) { cout<<"m["<<i<<"]="<<m[i]<<endl; }
// Формулы (7) for(int i=1;i<k;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["<<i-1<<"](x)= ("<< S[i-1][3] <<")*(x - (" <<p[i-1][0]<<"))^3 " <<" + ("<< S[i-1][2] <<")*(x - (" <<p[i-1][0]<<"))^2 + ("<< S[i-1][1] <<")*(x - ("<<p[i-1][0] <<")) + (" << S[i-1][0]<<")"<<endl; }
fstream outfile;
outfile.open("spline.txt", ios::out); //ios::app ios::in
outfile<<k<<" \n"; for(int i=0;i<k;i++) { outfile<<(p[i][0])<<" "; } outfile<<"\n"; for(int i=1;i<k;i++) { outfile<<S[i-1][3]<<" "<<S[i-1][2]<<" "<<S[i-1][1]<<" "<<S[i-1][0]<<" \n"; }
outfile.close();
system("pause"); return 0; } |
Двумерный кубический сплайн
|
#include <iostream> #include <cmath> #include <fstream>
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<nn;++i) { z[i]=new double[kk]; }
file.open("in1.txt"); for (int i=0; i < nn; ++i) { for(int j=0;j<kk;++j) { file>>z[i][j]; } } file.close ();
for (int i=0; i < nn; ++i) { for(int j=0;j<kk;++j) { cout<<z[i][j]<<" "; } cout<<endl; }
double** p= new double*[nn]; for(int i=0;i<nn;++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<nn;++i) { M[i]=new double[nn]; }
double** S= new double*[nn]; for(int i=0;i<nn;++i) { S[i]=new double[4]; }
//cout<<nn<<endl; //cout<<kk<<endl;
// 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"<<endl; cin>>x0>>xn; cin>>y0>>yn;
cout<<"Enter initial conditions for x , S'("<<x0<<") and S'("<<xn<<")"<<endl; cin>>sx0>>sxn;
cout<<"Enter initial conditions for y , S'("<<y0<<") and S'("<<yn<<")"<<endl; cin>>sy0>>syn;
///////////////// x
hx=(xn-x0)/(nn-1); //cout<<hx<<endl;
x[0]=x0; double temp=x0; for(int i=1;i<nn;i++) { temp=temp+hx; x[i]=temp; }
for(int imt=0;imt<kk;imt++) {
for(int i=0;i<nn;i++) { p[i][0]=x[i]; p[i][1]=z[i][imt]; // cout<<"("<<p[i][0]<<","<<p[i][1]<<")"<<endl; }
for(int i=1;i<nn;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<<h[i-1]<<endl; //cout<<d[i-1]<<endl; }
for(int i=1;i<nn-1;i++) { u[i]=6*(d[i]-d[i-1]); // формула (3) //cout<<u[i]<<endl; } // N=nn-1; // формула (4) if(nn==4) { for(int i=1;i<N;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<N;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<nn;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<nn;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<nn;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<<n<<endl;
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<nn;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<<x_point[i]<<endl; }
mm=0; for(int i=1;i<nn;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<<zx[mm][0]<<endl; mm=mm+1; } } zx[mm][imt]=z[nn-1][imt]; }
/////////////////// double** py= new double*[kk]; for(int i=0;i<kk;++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<kk;++i) { M1[i]=new double[kk]; }
double** S1= new double*[kk]; for(int i=0;i<kk;++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<<hy<<endl;
y[0]=y0; temp=y0;
for(int i=1;i<kk;i++) { temp=temp+hy; y[i]=temp; //cout<<y[i]<<endl; }
for(int jmt=0;jmt<=jj;jmt++) {
for(int i=0;i<kk;i++) { py[i][0]=y[i]; py[i][1]=zx[jmt][i]; //jmt // cout<<"("<<py[i][0]<<","<<py[i][1]<<")"<<endl; }
for(int i=1;i<kk;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<<h1[i-1]<<endl; //cout<<d1[i-1]<<endl; }
for(int i=1;i<kk-1;i++) { u1[i]=6*(d1[i]-d1[i-1]); //cout<<u1[i]<<endl; }
N=kk-1; // if(kk==4) { for(int i=1;i<N;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<N;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<kk;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<kk;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<kk;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<<n<<endl;
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<kk;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<kk;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<<x_point[i]<<" "; } outfile<<"\n"; for(int i=0;i<=jy;i++) { outfile<<y_point[i]<<" "; } outfile<<"\n";
for(int i=0;i<=jj;i++) { for(int j=0;j<=jy;j++) { outfile<<zxy[i][j]<<" "; } outfile<<"\n"; }
outfile<< nn <<" "<< kk <<" \n";
for(int i=0;i<nn;i++) { outfile<<x[i]<<" "; } outfile<<"\n"; for(int i=0;i<kk;i++) { outfile<<y[i]<<" "; } outfile<<"\n";
for(int i=0;i<nn;i++) { for(int j=0;j<kk;j++) { outfile<<z[i][j]<<" "; } outfile<<"\n"; }
outfile.close(); cout<<"ok"<<endl;
///////////// delete []d; delete []m; delete []u; delete []h; delete []x;
for(int i=0;i<nn;i++) delete [] z[i]; delete [] z;
for(int i=0;i<nn;i++) delete [] S[i]; delete [] S;
for(int i=0;i<nn;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 Кубический сплайн в формулах