利用c++编程三次样条插值,题目如下
发布网友
发布时间：2022-05-06 18:54
共1个回答
热心网友
时间：2023-09-11 06:39
#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;i<n;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;i<n-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;i<360*R/a;i++)
{
for(j=0;j<a/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;i<PI/2/a;i++)
{
for(j=0;j<a/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;
}
