三弯矩法

/*=======================================================================*/
#include <stdio.h>
#include <math.h>
///////////////////////////////////////////////////////////////////////////////#43;1]=B[k+1]-A[k+1]*C[k];
D[k+1]=D[k+1]-A[k+1]*D[k];
}
if(fabs(B[m-1])+1.0==1.0)
printf("fail at the second step\n");
D[m-1]=D[m-1]/B[m-1];
/*for(u=0;u<=m-1;u++)
printf("D[%d] is %f\n",u,D[u]);
for(u=0;u<=m-1;u++)
printf("B[%d] is %f\n",u,B[u]);*/
for(j=0;j<=m-1;j++)
{
printf("A[%d] is %f, B[%d] is %f, C[%d] is %f, D[%d] is %f\n ",j,A[j],j,B[j],j,C[j],j,D[j]);
}
for(i=0; i<=n;i++) printf("d[%d] is %f\n",i,d[i]);
printf("finish the calculation\n");
/* 追赶法求解M矩阵 */
m=n-1;
A[0]=0;
//C[0]=1;
double b3[NUM],b1[NUM]; //最终计算式系数
double xx=5.6; //输入节点
double Sxx=0.000000001; //插值结果
//分段函数的形式为 Si(x) = a3[i] * {x(i+1) - x}^3 + a1[i] * {x(i+1) - x} +
}
i=0;
while(xx>x[i+1]) i=i+1;
Sxx = a3[i] * (x[i+1] - xx)* (x[i+1] - xx)* (x[i+1] - xx) + a1[i] * (x[i+1] - xx) +b3[i] * (xx - x[i])* (xx - x[i])* (xx - x[i]) + b1[i] * (xx - x[i]);
for(k=m-2;k>=0;k--)
{D[k]=D[k]-C[k]*D[k+1];
}
for(i=0;i<=n;i++)
{M[i]=D[i];
printf("M[%d] is %f \n",i,M[i]);
}
for(i = 0;i<=n;i++)
{
{lambda[i]=h[i]/(h[i-1]+h[i]);
miu[i]=h[i-1]/(h[i-1]+h[i]);
d[i]=6*(f[i]-f[i-1])/(h[i-1]+h[i]);
// printf("miu[%d] is %f\n",i,miu[i]);
// printf("lambda[%d] is %f\n",i,lambda[i]);
printf("C[%d] is %f\n",u,C[u]);*/
for(k=0;k<=m-1;k++)
{B[k]=2;
D[k]=d[k+1];
}
D[0]=d[1]-miu[1]*ddy_begin;
D[m-1]=d[n-1]-lambda[n-1]*ddy_end;
/* 计算矩阵中间变量 */
for(i=0;i<=n-1;i++)
{h[i]=x[i+1]-x[i];
f[i]=(y[i+1]-y[i])/h[i];
// printf("h[%d] is %f\n",i,h[i]);
}
for(i=1;i<=n-1;i++)
printf("xx is %f\n",xx);
printf("result is %f\n",Sxx);
printf("end\n\n");
//fclose(fp);
_getch();
}
double B[NUM] = {0}; //追赶法中间变量
double C[NUM] = {0}; //追赶法中间量
double D[NUM] = {0}; //追赶法中间量
double a3[NUM],a1[NUM]; //最终计算式系数
double ddy_begin=0.0;
double ddy_end=0.0;
double dy[12];
double ddy[12];
void main()
{
double h[NUM] = {0}; //小区间的步长
double f[NUM] = {0}; //中间量
a3[i] = M[i] / (6 * h[i]);
a1[i] = (y[i] - M[i] * h[i] * h[i] / 6) / h[i];
b3[i] = M[i+1] / (6 * h[i]);
b1[i] = (y[i+1] - M[i+1] * h[i] * h[i] / 6) /h[i];
#define NUM 50 //定义样条数据区间个数最多为50个
//FILE *fp;
int n=11;
double x[12]={1,2,3,4,5,6,7,8,9,10,11,12};//{0.52, 8.0, 17.95, 28.65, 50.65, 104.6, 155.6, 260.7, 364.4, 468.0, 507.0, 520.0};//x[4]={0.0, 1.0, 4.0, 5.0};
for(k=1;k<=m-1;k++)
A[k]=miu[k+1];
for(k=0;k<=m-2;k++)
C[k]=lambda[k+1];
/*for(u=1;u<=m-1;u++)
printf("A[%d] is %f\n",u,A[u]);
for(u=0;u<=m-2;u++)
double y[12]={2,4,6,8,10,12,14,16,18,20,22,24};//{5.28794, 13.8400, 20.2000, 24.9000, 31.1000, 35.5000, 35.6000, 31.000, 20.9000, 7.80000, 1.50000, 0.20000};//y[4]={0.0, -2.0, -8.0, -4.0};
double miu[NUM] = {0}; //中间量
double lambda[NUM] = {0}; //中间量
double d[NUM] = {0}; //中间量
double M[NUM] = {0}; //M矩阵
double A[NUM] = {0}; //追赶法中间量
// b3[i] * {x - x(i)}^3 + b1[i] * {x - x(i)}
//xi为x[i]的值，xi_1为x[i+1]的值
int i = 0;
int m,k,u,j; //追赶法计数变量
printf("sv___111___\n");
}
printf("finish the initialization\n");
for(k=0;k<=m-2;k++)
{
if(fabs(B[k])+1.0==1.0)
printf("fail at the first step\n");
C[k]=C[k]/B[k];