建站
咨詢
售后
建站咨詢
新聞中心
求助：求三次樣條插值函數(shù)的C++程序
#includeiostream
成都創(chuàng)新互聯(lián)公司從2013年成立，是專業(yè)互聯(lián)網(wǎng)技術服務公司，擁有項目成都網(wǎng)站制作、做網(wǎng)站網(wǎng)站策劃，項目實施與項目整合能力。我們以讓每一個夢想脫穎而出為使命，1280元八公山做網(wǎng)站,已為上家服務,為八公山各地企業(yè)和個人服務,聯(lián)系電話:028-86922220
#includeiomanip
using
namespace
std;
const
int
MAX
=
50;
float
x[MAX],
y[MAX],
h[MAX];
float
c[MAX],
a[MAX],
fxym[MAX];
float
f(int
x1,
int
x2,
int
x3){
float
a
=
(y[x3]
-
y[x2])
/
(x[x3]
-
x[x2]);
float
b
=
(y[x2]
-
y[x1])
/
(x[x2]
-
x[x1]);
return
(a
-
b)/(x[x3]
-
x[x1]);
}
//求差分
void
cal_m(int
n){
//用追趕法求解出彎矩向量M……
float
B[MAX];
B[0]
=
c[0]
/
2;
for(int
i
=
1;
i
n;
i++)
B[i]
=
c[i]
/
(2
-
a[i]*B[i-1]);
fxym[0]
=
fxym[0]
/
2;
for(i
=
1;
i
=
n;
i++)
fxym[i]
=
(fxym[i]
-
a[i]*fxym[i-1])
/
(2
-
a[i]*B[i-1]);
for(i
=
n-1;
i
=
0;
i--)
fxym[i]
=
fxym[i]
-
B[i]*fxym[i+1];
}
void
printout(int
n);
int
main(){
int
n,i;
char
ch;
do{
cout"Please
put
in
the
number
of
the
dots:";
cinn;
for(i
=
0;
i
=
n;
i++){
cout"Please
put
in
X"i':';
cinx[i];
//coutendl;
cout"Please
put
in
Y"i':';
ciny[i];
//coutendl;
}
for(i
=
0;
i
n;
i++)
//求
步長
h[i]
=
x[i+1]
-
x[i];
cout"Please
輸入邊界條件\n
1:
已知兩端的一階導數(shù)\n
2:兩端的二階導數(shù)已知\n
默認:自然邊界條件\n";
int
t;
float
f0,
f1;
cint;
switch(t){
case
1:cout"Please
put
in
Y0\'
Y"n"\'\n";
cinf0f1;
c[0]
=
1;
a[n]
=
1;
fxym[0]
=
6*((y[1]
-
y[0])
/
(x[1]
-
x[0])
-
f0)
/
h[0];
fxym[n]
=
6*(f1
-
(y[n]
-
y[n-1])
/
(x[n]
-
x[n-1]))
/
h[n-1];
break;
case
2:cout"Please
put
in
Y0\"
Y"n"\"\n";
cinf0f1;
c[0]
=
a[n]
=
0;
fxym[0]
=
2*f0;
fxym[n]
=
2*f1;
break;
default:cout"不可用\n";//待定
};//switch
for(i
=
1;
i
n;
i++)
fxym[i]
=
6
*
f(i-1,
i,
i+1);
for(i
=
1;
i
n;
i++){
a[i]
=
h[i-1]
/
(h[i]
+
h[i-1]);
c[i]
=
1
-
a[i];
}
a[n]
=
h[n-1]
/
(h[n-1]
+
h[n]);
cal_m(n);
cout"\n輸出三次樣條插值函數(shù)：\n";
printout(n);
cout"Do
you
to
have
anther
try
?
y/n
:";
cinch;
}while(ch
==
'y'
||
ch
==
'Y');
return
0;
}
void
printout(int
n){
coutsetprecision(6);
for(int
i
=
0;
i
n;
i++){
couti+1":
["x[i]"
,
"x[i+1]"]\n""\t";
/*
coutfxym[i]/(6*h[i])"
*
("x[i+1]"
-
x)^3
+
""
*
(x
-
"x[i]")^3
+
"
(y[i]
-
fxym[i]*h[i]*h[i]/6)/h[i]"
*
("x[i+1]"
-
x)
+
"
(y[i+1]
-
fxym[i+1]*h[i]*h[i]/6)/h[i]"(x
-
"x[i]")\n";
coutendl;*/
float
t
=
fxym[i]/(6*h[i]);
if(t
0)coutt"*("x[i+1]"
-
x)^3";
else
cout-t"*(x
-
"x[i+1]")^3";
t
=
fxym[i+1]/(6*h[i]);
if(t
0)cout"
+
"t"*(x
-
"x[i]")^3";
else
cout"
-
"-t"*(x
-
"x[i]")^3";
cout"\n\t";
t
=
(y[i]
-
fxym[i]*h[i]*h[i]/6)/h[i];
if(t
0)cout"+
"t"*("x[i+1]"
-
x)";
else
cout"-
"-t"*("x[i+1]"
-
x)";
t
=
(y[i+1]
-
fxym[i+1]*h[i]*h[i]/6)/h[i];
if(t
0)cout"
+
"t"*(x
-
"x[i]")";
else
cout"
-
"-t"*(x
-
"x[i]")";
coutendlendl;
}
coutendl;
}
三次樣條插值用c語言具體怎么做
void SPL(int n, double *x, double *y, int ni, double *xi, double *yi)； 是你所要。
已知 n 個點 x,y; x 必須已按順序排好。要插值 ni 點，橫坐標 xi[], 輸出 yi[]。
程序里用double 型，保證計算精度。
SPL調用現(xiàn)成的程序。
現(xiàn)成的程序很多。端點處理方法不同，結果會有不同。想同matlab比較，你需 嘗試 調用 spline（）函數(shù) 時，令 end1 為 1, 設 slope1 的值，令 end2 為 1 設 slope2 的值。
誰來幫幫我啊:運用三次樣條插值函數(shù)法對一條已知的折線進行圓滑（編程實現(xiàn)）。
matlab算么？
function output = spline(x,y,xx)
%SPLINE Cubic spline data interpolation.
% PP = SPLINE(X,Y) provides the piecewise polynomial form of the
% cubic spline interpolant to the data values Y at the data sites X,
% for use with the evaluator PPVAL and the spline utility UNMKPP.
% X must be a vector.
% If Y is a vector, then Y(j) is taken as the value to be matched at X(j),
% hence Y must be of the same length as X -- see below for an exception
% to this.
% If Y is a matrix or ND array, then Y(:,...,:,j) is taken as the value to
% be matched at X(j), hence the last dimension of Y must equal length(X) --
% see below for an exception to this.
%
% YY = SPLINE(X,Y,XX) is the same as YY = PPVAL(SPLINE(X,Y),XX), thus
% providing, in YY, the values of the interpolant at XX. For information
% regarding the size of YY see PPVAL.
%
% Ordinarily, the not-a-knot end conditions are used. However, if Y contains
% two more values than X has entries, then the first and last value in Y are
% used as the endslopes for the cubic spline. If Y is a vector, this
% means:
% f(X) = Y(2:end-1), Df(min(X))=Y(1), Df(max(X))=Y(end).
% If Y is a matrix or N-D array with SIZE(Y,N) equal to LENGTH(X)+2, then
% f(X(j)) matches the value Y(:,...,:,j+1) for j=1:LENGTH(X), then
% Df(min(X)) matches Y(:,:,...:,1) and Df(max(X)) matches Y(:,:,...:,end).
%
% Example:
% This generates a sine-like spline curve and samples it over a finer mesh:
% x = 0:10; y = sin(x);
% xx = 0:.25:10;
% yy = spline(x,y,xx);
% plot(x,y,'o',xx,yy)
%
% Example:
% This illustrates the use of clamped or complete spline interpolation where
% end slopes are prescribed. In this example, zero slopes at the ends of an
% interpolant to the values of a certain distribution are enforced:
% x = -4:4; y = [0 .15 1.12 2.36 2.36 1.46 .49 .06 0];
% cs = spline(x,[0 y 0]);
% xx = linspace(-4,4,101);
% plot(x,y,'o',xx,ppval(cs,xx),'-');
%
% Class support for inputs x, y, xx:
% float: double, single
%
% See also INTERP1, PPVAL, UNMKPP, MKPP, SPLINES (The Spline Toolbox).
% Carl de Boor 7-2-86
% Copyright 1984-2004 The MathWorks, Inc.
% $Revision: 5.18.4.3 $ $Date: 2004/03/02 21:48:06 $
output=[];
% Check that data are acceptable and, if not, try to adjust them appropriately
[x,y,sizey,endslopes] = chckxy(x,y);
n = length(x); yd = prod(sizey);
% Generate the cubic spline interpolant in ppform
dd = ones(yd,1); dx = diff(x); divdif = diff(y,[],2)./dx(dd,:);
if n==2
if isempty(endslopes) % the interpolant is a straight line
pp=mkpp(x,[divdif y(:,1)],sizey);
else % the interpolant is the cubic Hermite polynomial
pp = pwch(x,y,endslopes,dx,divdif); pp.dim = sizey;
end
elseif n==3isempty(endslopes) % the interpolant is a parabola
y(:,2:3)=divdif;
y(:,3)=diff(divdif')'/(x(3)-x(1));
y(:,2)=y(:,2)-y(:,3)*dx(1);
pp = mkpp(x([1,3]),y(:,[3 2 1]),sizey);
else % set up the sparse, tridiagonal, linear system b = ?*c for the slopes
b=zeros(yd,n);
b(:,2:n-1)=3*(dx(dd,2:n-1).*divdif(:,1:n-2)+dx(dd,1:n-2).*divdif(:,2:n-1));
if isempty(endslopes)
x31=x(3)-x(1);xn=x(n)-x(n-2);
b(:,1)=((dx(1)+2*x31)*dx(2)*divdif(:,1)+dx(1)^2*divdif(:,2))/x31;
b(:,n)=...
(dx(n-1)^2*divdif(:,n-2)+(2*xn+dx(n-1))*dx(n-2)*divdif(:,n-1))/xn;
else
x31 = 0; xn = 0; b(:,[1 n]) = dx(dd,[2 n-2]).*endslopes;
end
dxt = dx(:);
c = spdiags([ [x31;dxt(1:n-2);0] ...
[dxt(2);2*[dxt(2:n-1)+dxt(1:n-2)];dxt(n-2)] ...
[0;dxt(2:n-1);xn] ],[-1 0 1],n,n);
% sparse linear equation solution for the slopes
mmdflag = spparms('autommd');
spparms('autommd',0);
s=b/c;
spparms('autommd',mmdflag);
% construct piecewise cubic Hermite interpolant
% to values and computed slopes
pp = pwch(x,y,s,dx,divdif); pp.dim = sizey;
end
if nargin==2, output = pp; else output = ppval(pp,xx); end
三次樣條插值C語言
我記得大三學的計算方法課上有，課后作業(yè)實現(xiàn)了的。不過在實驗室那個電腦上，如果你有條件的話先參考《數(shù)值分析》書上吧。
至于c語言和c++的區(qū)別，這個程序應該沒什么區(qū)別，反正都拿數(shù)組做。
(5)三次樣條和B樣條
首先什么是樣條？ 是 區(qū)間上的 個不同的點，當滿足如下兩個條件的時候， 就稱為一個 次的樣條函數(shù)
(1)
(2)
也就是說，在每個小區(qū)間上是 次多項式，并且整體是 階連續(xù)的。注意這里對樣條的定義并沒有要求在每個點 上函數(shù)值相等，如果真的滿足了 ，那么這就稱為樣條插值函數(shù)。
接下來有一個結論：
結論：對于 區(qū)間上由 個點所構建的所有的 次樣條函數(shù)張成的函數(shù)空間 ，其維數(shù) 。
這也就是說，對于任意一個樣條 ，它寫成基函數(shù)的形式應該是
???
而實際上，根據(jù)基函數(shù)的選擇不同，對應的樣條當然也不同，其中比較有名的是三次樣條和B樣條。
我們考慮 這一情況，并且是插值樣條，也就說 ，并且由于三次樣條要求二階連續(xù)，那么對于所有內部的節(jié)點 ，應該要求這些點處一階導數(shù)和二階導數(shù)應該左右相等。計算一下不難發(fā)現(xiàn)，這樣還缺少兩個條件，這里需要邊界條件，根據(jù)邊界條件的不同，插值樣條也不同。如自然樣條要求邊界點的二階導數(shù)為0。
構造三次樣條插值函數(shù)的方法如下，可以從每個區(qū)間的二階導數(shù)出發(fā)做一個線性插值，然后根據(jù)內部條件還有邊界條件構造方程組，最后解一個三對角的行列式。
B樣條的理論挺復雜的，在CAGD等領域是重點研究的方向。這里挖個坑，不寫了。
總結：插值先告一段落了。之前寫的都是一元函數(shù)的插值，實際上多元插值在研究領域中更重要，包括多元樣條，有限元等。從下一篇開始，寫數(shù)值積分。
三次樣條插值 C++程序
#includeiostream.h
#includeiomanip.h
#includemath.h
void main()
{
float a[37],b[37];
cout" ""度數(shù)"" ""sin(x)值"" ""一階導值"" ""二階導值"endl;
for(int i=0;i37;i++)
coutsetw(11)setprecision(3)
10*i
setw(11)setprecision(3)
sin(i*31.4/180)
setw(11)setprecision(3)
cos(i*31.4/180)
setw(11)setprecision(3)
-sin(i*31.4/180)
endl;
cout"一個周期內的積分值：0"endl;
}
名稱欄目：三次樣條函數(shù)c語言,c語言三次樣條插值函數(shù)
文章來源：http://www.dlmjj.cn/article/hojphh.html
