13518219792
#includeiostream
成都创新互联公司从2013年成立，是专业互联网技术服务公司，拥有项目成都网站制作、做网站网站策划，项目实施与项目整合能力。我们以让每一个梦想脱颖而出为使命，1280元八公山做网站,已为上家服务,为八公山各地企业和个人服务,联系电话: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:
已知两端的一阶导数\n
2:两端的二阶导数已知\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输出三次样条插值函数：\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;
}
void SPL(int n, double *x, double *y, int ni, double *xi, double *yi)； 是你所要。
已知 n 个点 x,y; x 必须已按顺序排好。要插值 ni 点，横坐标 xi[], 输出 yi[]。
程序里用double 型，保证计算精度。
SPL调用现成的程序。
现成的程序很多。端点处理方法不同，结果会有不同。想同matlab比较，你需 尝试 调用 spline（）函数 时，令 end1 为 1, 设 slope1 的值，令 end2 为 1 设 slope2 的值。
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语言和c++的区别，这个程序应该没什么区别，反正都拿数组做。
首先什么是样条？ 是 区间上的 个不同的点，当满足如下两个条件的时候， 就称为一个 次的样条函数
(1)
(2)
也就是说，在每个小区间上是 次多项式，并且整体是 阶连续的。注意这里对样条的定义并没有要求在每个点 上函数值相等，如果真的满足了 ，那么这就称为样条插值函数。
接下来有一个结论：
结论：对于 区间上由 个点所构建的所有的 次样条函数张成的函数空间 ，其维数 。
这也就是说，对于任意一个样条 ，它写成基函数的形式应该是
而实际上，根据基函数的选择不同，对应的样条当然也不同，其中比较有名的是三次样条和B样条。
我们考虑 这一情况，并且是插值样条，也就说 ，并且由于三次样条要求二阶连续，那么对于所有内部的节点 ，应该要求这些点处一阶导数和二阶导数应该左右相等。计算一下不难发现，这样还缺少两个条件，这里需要边界条件，根据边界条件的不同，插值样条也不同。如自然样条要求边界点的二阶导数为0。
构造三次样条插值函数的方法如下，可以从每个区间的二阶导数出发做一个线性插值，然后根据内部条件还有边界条件构造方程组，最后解一个三对角的行列式。
B样条的理论挺复杂的，在CAGD等领域是重点研究的方向。这里挖个坑，不写了。
总结：插值先告一段落了。之前写的都是一元函数的插值，实际上多元插值在研究领域中更重要，包括多元样条，有限元等。从下一篇开始，写数值积分。
#includeiostream.h
#includeiomanip.h
#includemath.h
void main()
{
float a[37],b[37];
cout" ""度数"" ""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;
}
