三次样条插值用c语言具体怎么做

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 的值。

#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

}

#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=0i<=360/ai++)

{

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=0i<=360/ai++)

{

f1[i]=*(y1+i)

}

///打印///

printf("角度值\t|函数值\t|一阶导数值\t\n")

for(i=0i<=360/ai++)

{

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=1i<ni++)

{

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=0i<=ni++)

{

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=2i<n-1i++)

{

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-2i>=1i--)

{

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=0i<360*R/ai++)

{

for(j=0j<a/Rj++)

{

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=0i<PI/2/ai++)

{

for(j=0j<a/Rj++)

{

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

}

---