#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
class scyt
{private:
int i,n;
double s;//插入点的计算值
double X;//输入的插入点
double z1,z2;//起点和终点的一阶导数值
double *xx,*yy,*a,*b,*c,*f,*h;//b用来存主对角线上的值，
//a,c用来存下，上次对角线上的值
double *x,*y,*l,*u;//求解三对角方程所需元素
public:
void input();//数据输入
void qiudao();//求二阶导数
void sdj();//解三对角方程
void jscrd();//计算插入点的值
~scyt()//释放内存
{delete []xx,yy,a,b,c,h,x,y,l,u,f;}
};
void main()
{scyt myscyt;
myscyt.input();//读入x,y数据
myscyt.qiudao();//计算二阶导数
myscyt.jscrd();//计算插入点的值
}
void scyt::input()//读入x,y数据
{ifstream fin("ty.dat");
fin>>n;
xx=new double[n];
yy=new double[n];
for(i=0;i<n;i++)
{fin>>xx[i]>>yy[i];}
fin>>z1>>z2;
fin.close();
}
void scyt::qiudao()//求二阶导数
{a=new double[n-1];
b=new double[n];
c=new double[n-1];
f=new double[n];
h=new double[n-1];
for(i=0;i<n-1;i++)
{h[i]=xx[i+1]-xx[i];}
a[n-2]=1;
c[0]=1;
for(i=0;i<n-2;i++)
{a[i]=h[i]/(h[i]+h[i+1]);}
for(i=0;i<n;i++)
{b[i]=2;}
for(i=1;i<n-1;i++)
{c[i]=1-a[i-1];}
f[0]=(6/h[0])*((yy[1]-yy[0])-z1);
f[n-1]=(6/h[n-2])*(z2-(yy[n-1]-yy[n-2])/h[n-2]);
for(i=1;i<n-1;i++)
{f[i]=(6/(h[i-1]+h[i]))*((yy[i+1]-yy[i])/h[i]-\
(yy[i]-yy[i-1])/h[i-1]);}
sdj();
}
void scyt::sdj()//解三对角方程
{x=new double[n];
y=new double[n];
l=new double[n];
u=new double[n-1];
l[0]=b[0];
for(i=0;i<=n-2;i++)
{u[i]=c[i]/l[i];
l[i+1]=b[i+1]-a[i]*u[i];}//求下三角中的l[i]和上三角中的u[i]
y[0]=f[0]/l[0];
for(i=1;i<n;i++)
{y[i]=(f[i]-a[i-1]*y[i-1])/l[i];}
x[n-1]=y[n-1];
for(i=n-2;i>=0;i--)
{x[i]=y[i]-u[i]*x[i+1];}
cout<<n<<"个点的导数值分别是:"<<endl;
for(i=0;i<n;i++)
cout<<"x["<<i<<"]="<<x[i]<<" ";
}
void scyt::jscrd()
{cout<<"\n"<<"请输入插入点的值X,当X=1000时退出程序！"<<endl;
while(X!=1000)
{cout<<"X=:";
cin>>X;
if(X!=1000)
{for(i=0;i<n-1;i++)
{if(xx[i]<=X&&xx;[i+1]>=X)
break;
}
s=pow((xx[i+1]-X),3)*x[i]/(6*h[i])+pow\
((X-xx[i]),3)*x[i+1]/(6*h[i])+\
(xx[i+1]-X)*(yy[i]-h[i]*h[i]*x[i]/6)/h[i]\
+(X-xx[i])*(yy[i+1]-h[i]*h[i]*x[i+1]/6)/h[i];
cout<<"s["<<X<<"]="<<s<<endl;
}
}
}
0
106
2KB
2012-06-12
9
