#include stdio.h
创新互联建站是专业的乌审网站建设公司，乌审接单;提供网站建设、成都网站建设,网页设计,网站设计,建网站,PHP网站建设等专业做网站服务;采用PHP框架,可快速的进行乌审网站开发网页制作和功能扩展;专业做搜索引擎喜爱的网站,专业的做网站团队,希望更多企业前来合作!
#include math.h
double fun(double x,double x0,double x1,double x2,double y0,double y1,double y2)
{
double yx=0;
yx=y0*(x-x1)*(x-x2)/((x0-x1)*(x0-x2))+
y1*(x-x0)*(x-x2)/((x1-x0)*(x1-x2))+
y2*(x-x0)*(x-x1)/((x2-x0)*(x2-x1));//3点插值公式
return yx;
}
int main(int argc, char *argv[])
{
double x,x0,x1,x2,y0,y1,y2;
printf("输入待求值x：\n");
scanf("%lf",x);
x0=x-0.1;x1=x+0.1;x2=x+0.15;//需要输入3个插值点，即对应的x值和函数y值，这里简单计算的可以手动输入
y0=sin(x0);y1=sin(x1);y2=sin(x2);
printf("sin(%lf)=%lf-------fun(%lf)=%lf\n",x,sin(x),x,fun(x,x0,x1,x2,y0,y1,y2));
return 0;
}
我记得大三学的计算方法课上有，课后作业实现了的。不过在实验室那个电脑上，如果你有条件的话先参考《数值分析》书上吧。
至于c语言和c++的区别，这个程序应该没什么区别，反正都拿数组做。
三次样条插值在实际中有着广泛的应用，在计算机上也容易实现。下面介绍用计算机求取三样条插值函数S（x）的算法步骤：
（1）输入初始节点离散数据xi，yi（i＝0，1，…，n）；
（2）依据式（6－46），计算hi＝xi－xi－1，λi和Ri（i＝1，…，n－1）；
（3）根据实际问题，从式（6－49）、式（6－51）和式（6－53）中选择一类对应的边界条件，求取v0，w0，u0，R0，un，vn，wn，Rn；
（4）根据形成的方程组（6－54）的特点，选用追赶法、高斯法等解方程组，求出Mi（i＝0，1，2，…，n）；
（5）依据式（6－41）、式（6－42），计算插值点的三样条插值函数值和该点的导数值。
#includeiostream
#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;
}
#include"stdio.h"
#include"math.h"
const double PI = 3.141592654;
const double R=PI/180; //定义常量（角度转弧度的因子）。
double *m(double *y,int length,double j); ////m法////
double *jinsi(double y[],double g[],double a); ///求给定间隔的函数近似值///
double *jifen(double y[],double g[],double a,double b);//求第一象限内的积分
void main()
{
/////给定已知条件////
int a; //给定间隔（角度值a）
double y[361],*p,*q; //定义名为y的数组存储函数值，下标及对应角度
double* y1, * y2,*y3,*y4;
double f1[361],f2[361];
printf("请输入已知sin函数表的间隔a：");
scanf("%d",a); //就是h （等间距）
for(int i=0;i=360/a;i++)
{
if(i*a%180==0)
y[i]=0.0;
else y[i]=sin(R*i*a); //存储函数值
}
p=y;
y1=m(p,360/a+1,R*a); //用m法求一阶导(1~n-1)
for(i=0;i=360/a;i++)
{
f1[i]=*(y1+i);
}
///打印///
printf("角度值\t|函数值\t|一阶导数值\t\n");
for(i=0;i=360/a;i++)
{
if(i/10==0) printf("%d | %f | %f\n",i*a,y[i],f1[i]);
else if(i/100==0) printf("%d | %f | %f\n",i*a,y[i],f1[i]);
else printf("%d | %f | %f\n",i*a,y[i],f1[i]);
}
}
////m法////
double *m(double *y,int length,double j)
{
int n=length-1;
double g1[361];
for(int i=1;in;i++)
{
g1[i]=(*(y+i+1)-*(y+i-1))*3/(j*2);
}
g1[0]=1;
g1[n]=1;
g1[1]-=0.5*g1[0];
g1[n-1]-=0.5*g1[n];
double a[360];
double c[360];
double d[360];
for(i=0;i=n;i++)
{
a[i]=0.5;
d[i]=2;
c[i]=0.5;
}
c[1]=c[1]/d[1];
g1[1]=g1[1]/d[1];
double t;
for(i=2;in-1;i++)
{
t=d[i]-c[i-1]*a[i];
c[i]/=t;
g1[i]=(g1[i]-g1[i-1]*a[i])/t;
}
g1[n-1]=(g1[n-1]-g1[n-2]*a[n-1])/(d[n-1]-c[n-2]*a[n-1]);
for(i=n-2;i=1;i--)
{
g1[i]=g1[i]-c[i]*g1[i+1];
}
return g1;
}
///求近似函数值//
double *jinsi(double y[],double g[],double a)//a为弧度值
{
int i=0,j=0;
double f[361];
for(i=0;i360*R/a;i++)
{
for(j=0;ja/R;j++)
{
double x,s,a1,a2,b1,b2;
x=i*a+j*R;//变为弧度制
a1=(1+2*(x-a*i)/a)*(x-a*(i+1))*(x-a*(i+1))/a/a;
a2=(1-2*(x-a*(i+1))/a)*(x-a*i)*(x-a*i)/a/a;
b1=(x-a*i)*(x-a*(i+1))*(x-a*(i+1))/a/a;
b2=(x-a*(i+1))*(x-a*i)*(x-a*i)/a/a;
s=a1*y[i]+y[i+1]*a2+g[i]*b1*g[i+1]*b2; //间隔为10的样条表达式
f[i*10+j]=s;
}
}
f[360]=0;
return f;
}
double *jifen(double y[],double g[],double a,double b)//a为弧度值
{ int i=0,j=0;
double sum=0;
for(i=0;iPI/2/a;i++)
{
for(j=0;ja/R;j++)
{
double x,s,a1,a2,b1,b2;
x=i*a+j*R;//变为弧度制
a1=(1+2*(x-a*i)/a)*(x-a*(i+1))*(x-a*(i+1))/a/a;
a2=(1-2*(x-a*(i+1))/a)*(x-a*i)*(x-a*i)/a/a;
b1=(x-a*i)*(x-a*(i+1))*(x-a*(i+1))/a/a;
b2=(x-a*(i+1))*(x-a*i)*(x-a*i)/a/a;
s=a1*y[i]+y[i+1]*a2+g[i]*b1*g[i+1]*b2; //间隔为10的样条表达式
sum+=s*R;
}
}
printf("第一象限内的积分值为：%f/n",sum);
return 0;
}
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 的值。
重庆分公司
重庆分公司
