
x=0:10
y=sin(x)
xx=0:.25:10
yy=spline(x,y,xx)
plot(x,y,'o',xx,yy)
另外fnpltcsapi这两个函数也是三次样条插值函数,具体你可以help一下!
现在电脑上没有matlab,一会给你程序,呵呵!
#include<iostream>#include<iomanip>
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:"
cin>>n
for(i
=
0
i
<=
n
i++){
cout<<"Please
put
in
X"<<i<<':'
cin>>x[i]
//cout<<endl
cout<<"Please
put
in
Y"<<i<<':'
cin>>y[i]
//cout<<endl
}
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
cin>>t
switch(t){
case
1:cout<<"Please
put
in
Y0\'
Y"<<n<<"\'\n"
cin>>f0>>f1
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"
cin>>f0>>f1
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
:"
cin>>ch
}while(ch
==
'y'
||
ch
==
'Y')
return
0
}
void
printout(int
n){
cout<<setprecision(6)
for(int
i
=
0
i
<
n
i++){
cout<<i+1<<":
["<<x[i]<<"
,
"<<x[i+1]<<"]\n"<<"\t"
/*
cout<<fxym[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"
cout<<endl*/
float
t
=
fxym[i]/(6*h[i])
if(t
>
0)cout<<t<<"*("<<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]<<")"
cout<<endl<<endl
}
cout<<endl
}
欢迎分享,转载请注明来源:内存溢出
微信扫一扫
支付宝扫一扫
评论列表(0条)