3,881
社区成员
发帖
与我相关
我的任务
分享求空间插值C++代码！
有离散点数据和坐标，插值成可以画等值线的格点文件，保存在txt文件中，只要求插值部分，不用画图，可以生成格点文件即可，谢谢！！克里金插值什么的都可以！
...全文
请发表友善的回复…
发表回复
mujiok2003 2013-08-18
- 打赏
- 举报
这叫样条差值[/quote]
是把二维离散点网格化吗[/quote]
自己也看看。
mujiok2003 2013-08-18
- 打赏
- 举报
这叫样条差值
derekrose 2013-08-18
- 打赏
- 举报
网上应该一大堆吧
yaoyunz 2013-08-18
- 打赏
- 举报
但是每一步都是做什么的呢 根本看不懂啊？[/quote]
读懂代码需要找到借助“数值计算”和高等数学。 完整代码
#include <iostream>
using namespace std;
//from the book "Numeric recipes"
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
int i,k;
double p,qn,sig,un,*u;
u = (double *) malloc((size_t) (n*sizeof(double)));
if(yp1 > 0.99e30)
{
// "Natural" 1st derivative...
y2[0]=u[0]=0.0;
}else
{
// Specified 1st derivative...
y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0])/(x[1]-x[0]) - yp1);
}
for(i=1;i<=n-2;i++)
{
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
p=sig*y2[i-1] + 2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
}
if(ypn > 0.99e30)
{
qn = un = 0.0;
}else
{
qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2]+1.0);
for(k=n-2;k>=0;k--)
{
y2[k] = y2[k] * y2[k+1] + u[k];
}
free(u);
}
//from the book "Numeric recipes"
double splint(double xa[], double ya[], double y2a[], int n, double x)
{
int klo,khi,k;
double h,b,a,y;
klo = 0;
khi = n-1;
while((khi-klo) > 1)
{
k=(khi+klo) >> 1;
if(xa[k] > x)
{
khi = k;
}else
{
klo = k;
}
}
h = xa[khi] - xa[klo];
// if(h == 0.0) ERROR !!!
a = (xa[khi]-x)/h;
b = (x - xa[klo])/h;
y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return (y);
}
#define N_POINTS 10
int main()
{
int n;
double yp1, ypn, x[N_POINTS], y[N_POINTS], y2[N_POINTS];
double input, output;
int i;
n = N_POINTS;
//Natural cubic spline
yp1 = 1.0e30;
ypn = 1.0e30;
//y = x*x
for(i=0;i<N_POINTS;++i)
{
x[i] = i;
y[i] = i*i;
}
spline(x, y, n, yp1, ypn, y2);
for(i=0;i<N_POINTS;++i)
{
input = (double) i;
output = splint(x, y, y2, n, input);
cout << "input output " << input << " " << output << "\n";
if(i != (N_POINTS - 1))
{
input = (double) i + 0.5;
output = splint(x, y, y2, n, input);
cout << "input output " << input << " " << output << "\n";
}
}
return 0;
}
mujiok2003 2013-08-16
- 打赏
- 举报
但是每一步都是做什么的呢 根本看不懂啊？[/quote]
读懂代码需要找到借助“数值计算”和高等数学。 完整代码
#include <iostream>
using namespace std;
//from the book "Numeric recipes"
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
int i,k;
double p,qn,sig,un,*u;
u = (double *) malloc((size_t) (n*sizeof(double)));
if(yp1 > 0.99e30)
{
// "Natural" 1st derivative...
y2[0]=u[0]=0.0;
}else
{
// Specified 1st derivative...
y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0])/(x[1]-x[0]) - yp1);
}
for(i=1;i<=n-2;i++)
{
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
p=sig*y2[i-1] + 2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
}
if(ypn > 0.99e30)
{
qn = un = 0.0;
}else
{
qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2]+1.0);
for(k=n-2;k>=0;k--)
{
y2[k] = y2[k] * y2[k+1] + u[k];
}
free(u);
}
//from the book "Numeric recipes"
double splint(double xa[], double ya[], double y2a[], int n, double x)
{
int klo,khi,k;
double h,b,a,y;
klo = 0;
khi = n-1;
while((khi-klo) > 1)
{
k=(khi+klo) >> 1;
if(xa[k] > x)
{
khi = k;
}else
{
klo = k;
}
}
h = xa[khi] - xa[klo];
// if(h == 0.0) ERROR !!!
a = (xa[khi]-x)/h;
b = (x - xa[klo])/h;
y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return (y);
}
#define N_POINTS 10
int main()
{
int n;
double yp1, ypn, x[N_POINTS], y[N_POINTS], y2[N_POINTS];
double input, output;
int i;
n = N_POINTS;
//Natural cubic spline
yp1 = 1.0e30;
ypn = 1.0e30;
//y = x*x
for(i=0;i<N_POINTS;++i)
{
x[i] = i;
y[i] = i*i;
}
spline(x, y, n, yp1, ypn, y2);
for(i=0;i<N_POINTS;++i)
{
input = (double) i;
output = splint(x, y, y2, n, input);
cout << "input output " << input << " " << output << "\n";
if(i != (N_POINTS - 1))
{
input = (double) i + 0.5;
output = splint(x, y, y2, n, input);
cout << "input output " << input << " " << output << "\n";
}
}
return 0;
}
yaoyunz 2013-08-16
- 打赏
- 举报
我一般用“数值计算”那本书的代码
//from the book "Numeric recipes"
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
int i,k;
double p,qn,sig,un,*u;
u = (double *) malloc((size_t) (n*sizeof(double)));
if(yp1 > 0.99e30)
{
// "Natural" 1st derivative...
y2[0]=u[0]=0.0;
}else
{
// Specified 1st derivative...
y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0])/(x[1]-x[0]) - yp1);
}
for(i=1;i<=n-2;i++)
{
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
p=sig*y2[i-1] + 2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
}
if(ypn > 0.99e30)
{
qn = un = 0.0;
}else
{
qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2]+1.0);
for(k=n-2;k>=0;k--)
{
y2[k] = y2[k] * y2[k+1] + u[k];
}
free(u);
}
//from the book "Numeric recipes"
double splint(double xa[], double ya[], double y2a[], int n, double x)
{
int klo,khi,k;
double h,b,a,y;
klo = 0;
khi = n-1;
while((khi-klo) > 1)
{
k=(khi+klo) >> 1;
if(xa[k] > x)
{
khi = k;
}else
{
klo = k;
}
}
h = xa[khi] - xa[klo];
// if(h == 0.0) ERROR !!!
a = (xa[khi]-x)/h;
b = (x - xa[klo])/h;
y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return (y);
}
mujiok2003 2013-08-15
- 打赏
- 举报
我一般用“数值计算”那本书的代码
//from the book "Numeric recipes"
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
int i,k;
double p,qn,sig,un,*u;
u = (double *) malloc((size_t) (n*sizeof(double)));
if(yp1 > 0.99e30)
{
// "Natural" 1st derivative...
y2[0]=u[0]=0.0;
}else
{
// Specified 1st derivative...
y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0])/(x[1]-x[0]) - yp1);
}
for(i=1;i<=n-2;i++)
{
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
p=sig*y2[i-1] + 2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
}
if(ypn > 0.99e30)
{
qn = un = 0.0;
}else
{
qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2]+1.0);
for(k=n-2;k>=0;k--)
{
y2[k] = y2[k] * y2[k+1] + u[k];
}
free(u);
}
//from the book "Numeric recipes"
double splint(double xa[], double ya[], double y2a[], int n, double x)
{
int klo,khi,k;
double h,b,a,y;
klo = 0;
khi = n-1;
while((khi-klo) > 1)
{
k=(khi+klo) >> 1;
if(xa[k] > x)
{
khi = k;
}else
{
klo = k;
}
}
h = xa[khi] - xa[klo];
// if(h == 0.0) ERROR !!!
a = (xa[khi]-x)/h;
b = (x - xa[klo])/h;
y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return (y);
}
yaoyunz 2013-08-15
- 打赏
- 举报
我一般用“数值计算”那本书的代码
//from the book "Numeric recipes"
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
int i,k;
double p,qn,sig,un,*u;
u = (double *) malloc((size_t) (n*sizeof(double)));
if(yp1 > 0.99e30)
{
// "Natural" 1st derivative...
y2[0]=u[0]=0.0;
}else
{
// Specified 1st derivative...
y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0])/(x[1]-x[0]) - yp1);
}
for(i=1;i<=n-2;i++)
{
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
p=sig*y2[i-1] + 2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
}
if(ypn > 0.99e30)
{
qn = un = 0.0;
}else
{
qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2]+1.0);
for(k=n-2;k>=0;k--)
{
y2[k] = y2[k] * y2[k+1] + u[k];
}
free(u);
}
//from the book "Numeric recipes"
double splint(double xa[], double ya[], double y2a[], int n, double x)
{
int klo,khi,k;
double h,b,a,y;
klo = 0;
khi = n-1;
while((khi-klo) > 1)
{
k=(khi+klo) >> 1;
if(xa[k] > x)
{
khi = k;
}else
{
klo = k;
}
}
h = xa[khi] - xa[klo];
// if(h == 0.0) ERROR !!!
a = (xa[khi]-x)/h;
b = (x - xa[klo])/h;
y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return (y);
}
mujiok2003 2013-08-14
- 打赏
- 举报
我一般用“数值计算”那本书的代码
//from the book "Numeric recipes"
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
int i,k;
double p,qn,sig,un,*u;
u = (double *) malloc((size_t) (n*sizeof(double)));
if(yp1 > 0.99e30)
{
// "Natural" 1st derivative...
y2[0]=u[0]=0.0;
}else
{
// Specified 1st derivative...
y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0])/(x[1]-x[0]) - yp1);
}
for(i=1;i<=n-2;i++)
{
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
p=sig*y2[i-1] + 2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
}
if(ypn > 0.99e30)
{
qn = un = 0.0;
}else
{
qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2]+1.0);
for(k=n-2;k>=0;k--)
{
y2[k] = y2[k] * y2[k+1] + u[k];
}
free(u);
}
//from the book "Numeric recipes"
double splint(double xa[], double ya[], double y2a[], int n, double x)
{
int klo,khi,k;
double h,b,a,y;
klo = 0;
khi = n-1;
while((khi-klo) > 1)
{
k=(khi+klo) >> 1;
if(xa[k] > x)
{
khi = k;
}else
{
klo = k;
}
}
h = xa[khi] - xa[klo];
// if(h == 0.0) ERROR !!!
a = (xa[khi]-x)/h;
b = (x - xa[klo])/h;
y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return (y);
}
mujiok2003 2013-08-14
- 打赏
- 举报
我一般用“数值计算”那本书的代码
//from the book "Numeric recipes"
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
int i,k;
double p,qn,sig,un,*u;
u = (double *) malloc((size_t) (n*sizeof(double)));
if(yp1 > 0.99e30)
{
// "Natural" 1st derivative...
y2[0]=u[0]=0.0;
}else
{
// Specified 1st derivative...
y2[0] = -0.5;
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0])/(x[1]-x[0]) - yp1);
}
for(i=1;i<=n-2;i++)
{
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
p=sig*y2[i-1] + 2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i] = (6.0*u[i]/(x[i+1]-x[i-1]) - sig*u[i-1])/p;
}
if(ypn > 0.99e30)
{
qn = un = 0.0;
}else
{
qn = 0.5;
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
}
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2]+1.0);
for(k=n-2;k>=0;k--)
{
y2[k] = y2[k] * y2[k+1] + u[k];
}
free(u);
}
//from the book "Numeric recipes"
double splint(double xa[], double ya[], double y2a[], int n, double x)
{
int klo,khi,k;
double h,b,a,y;
klo = 0;
khi = n-1;
while((khi-klo) > 1)
{
k=(khi+klo) >> 1;
if(xa[k] > x)
{
khi = k;
}else
{
klo = k;
}
}
h = xa[khi] - xa[klo];
// if(h == 0.0) ERROR !!!
a = (xa[khi]-x)/h;
b = (x - xa[klo])/h;
y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return (y);
}
yaoyunz 2013-08-14
- 打赏
- 举报
哪个能用啊
mujiok2003 2013-08-14
- 打赏
- 举报
