急求三次样条的算法C程序

湘军 · 发表于 2005-4-20 05:03:22

<

>帮我看一下，这程序那出错了，高手们帮我改一下，要是你们有更好的，请发出来，我急用！

<

>#include <stdio.h>
#include <math.h>
#include <memory.h> 
<

>void swap(double & a, double & b)
{
double temp;
temp=a;
a=b;
b=temp;
} 
bool GuassSolve(double *A,double *B,double *X,int n)
{
int i,j,k,row;
for(k=0;k<n-1;k++)
{
double max=0.0;
for(i=k;i<n;i++)
{
if(fabs(A[i*n+k])>max)
{
max=fabs(A[i*n+k]);
row=i;
}
}
if(max==0.0) return false;
if(row!=k)
{
for(j=0;j<n;j++)
swap(A[row*n+j],A[k*n+j]);
swap(B[row],B[k]);
}
for(i=k+1;i<n;i++)
{
double m=A[i*n+k]/A[k*n+k];
for(j=k+1;j<n;j++)
{
A[i*n+j]=A[i*n+j]-m*A[k*n+j];
}
B=B-m*B[k];
}
}
X[n-1]=B[n-1]/A[n*n-1];
for(k=n-2;k>=0;k--)
{
double sum=0.0;
for(j=k+1;j<n;j++)
sum+=A[k*n+j]*X[j];
X[k]=(B[k]-sum)/A[k*n+k];
}
return true;
} 
void Apend(double miu[],double rmd[],double d[],double x[],double y[],double h

[],
 double pf1,double pf2,
 double ppf1,double ppf2,
 int n,int flag)
{
switch(flag)
{
case 0:
miu[n-1]=1;
rmd[0]=1;
d[0]=6*((y[1]-y[0])/(x[1]-x[0])-pf1)/h[0];
d[n-1]=6*(pf2-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/h[n-2];
break;
case 1:
rmd[0]=miu[n-1]=0;
d[0]=2*ppf1;
d[n-1]=2*ppf2;
break;
case 2:
rmd[0]=rmd[n-1]=h[0]/(h[n-1]+h[0]);
miu[0]=miu[n-1]=1-rmd[n-1];
d[0]=d[n-1]=6*((y[1]-y[0])/(x[1]-x[0])-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/(h 
[0 
]+h[n-2]);
break;
}
} 
void T(double x[],double y[],int n)
{
int flag=0;
double pf1=3.0;
double pf2=-4.0;
double ppf1=0;
double ppf2=0;
double * h=new double [n];
double * miu=new double [n];
double * rmd=new double [n];
double * d=new double[n]; 
memset(d,0,n*sizeof(double));
memset(rmd,0,n*sizeof(double));
memset(miu,0,n*sizeof(double)); 
for(int i=1;i<n;i++)
{
h[i-1]=x-x[i-1];
} 
for(i=1;i<n-1;i++)
{
rmd=h/(h[i-1]+h);
miu=h[i-1]/(h[i-1]+h);
d=6*((y[i+1]-y)/(x[i+1]-x)-(y-y[i-1])/(x-x[i-1]))/(h[i-1]+h 
[i 
]);
} 
Apend(miu,rmd,d,x,y,h,pf1,pf2,ppf1,ppf2,n,0); 
double * A=new double [n*n],* M=new double [n]; 
memset(A,0,n*n*sizeof(double)); 
for(i=1;i<n-1;i++)
{
A[i*n+i]=2; 
A[i*n+i-1]=miu; 
A[i*n+i+1]=rmd;
} 
if(n==2)
{
A[0]=2;
A[n-2]=miu[0];
A[1]=rmd[0];
A[(n-1)*n+n-1]=2;
A[(n-1)*n+n-2]=miu[n-1];
A[(n-1)*n+1]=rmd[n-1];
}
else
{
A[0]=2;
A[1]=rmd[0];
A[(n-1)*n+n-1]=2;
A[(n-1)*n+n-2]=miu[n-1];
} 
GuassSolve(A,d,M,n); 
for(i=0;i<n;i++)
{
printf("M[%d]=%f\n",i+1,M);
} 
} 
void Initiate(double * x,double * y)
{
x[0]=27.7;
y[0]=4.1; 
x[1]=28;
y[1]=4.3; 
x[2]=29;
y[2]=4.1; 
x[3]=30;
y[3]=3.0;
} 
void main()
{
int n=4;
double * x=new double[n];
double * y=new double[n]; 
Initiate(x,y); 
T(x,y,n);
}

城市卖风 · 发表于 2005-4-21 17:43:26

<

>浩大一个工程啊？抱歉，无能为力，我有Ｃ＋＋的问题呢，不过有个贴子好象很多程序，你到１１页后看把！

sunjlee · 发表于 2005-5-3 16:07:35

一开始swap ()里 temp好象要指针型吧

		自动登录	找回密码
密码			注-册-帐-号