16.txt

追赶法解三对角线方程组

#include <math.h>
#include <malloc.h>
#include <stdio.h>
#include <stdlib.h>
#include <graphics.h>
int n,N;
/*追赶法解三对角线方程组*/
void RunAfterMethod(double *u,double *v,double *g,double *M){
   double *y,*b;int i;
   y=(double*)malloc(N*sizeof(double));
   b=(double*)malloc(N*sizeof(double));
   for(*b=*(b+n)=0,i=1;i<=n-1;i++)
      *(b+i)=*(v+i)/(2-*(u+i)**(b+i-1));
   for(*y=0,i=1;i<=n;i++)
      *(y+i)=(*(g+i)-*(u+i)**(y+i-1))/(2-*(u+i)**(b+i-1));
   for(*(M+n)=*M=0,i=n-1;i>0;i--)
      *(M+i)=*(y+i)-*(b+i)**(M+i+1);
   free(y),free(b);
}
/*自然样条插值函数*/
double S(double dx,double *x,double *y,double *h,double *M){
   int i;
   if(dx<*x)return(*y);
   if(dx>*(x+n))return(*(y+n));
   for(i=1;i<=n;i++)
      if(dx>=*(x+i-1)&&dx<=*(x+i))break;
   return *(M+i-1)*pow(*(x+i)-dx,3)/(6**(h+i))+
        *(M+i)*pow(dx-*(x+i-1),3)/(6**(h+i))+
        (*(y+i-1)-*(M+i-1)*pow(*(h+i),2)/6)*(*(x+i)-dx)/(*(h+i))+
        (*(y+i)-*(M+i)*pow(*(h+i),2)/6)*(dx-*(x+i-1))/(*(h+i));
}
int main(void){
   double *x,*y,*h,*u,*v,*g,*M,t1,t2,t=0;
   FILE *fi;
   int gd=DETECT,gm,i;
   fi=fopen("input.dat","r");
   fscanf(fi,"%d",&N);n=N-1;
   x=(double*)malloc(N*sizeof(double));
   y=(double*)malloc(N*sizeof(double));
   h=(double*)malloc(N*sizeof(double));
   u=(double*)malloc(N*sizeof(double));
   v=(double*)malloc(N*sizeof(double));
   g=(double*)malloc(N*sizeof(double));
   M=(double*)malloc(N*sizeof(double));
   for(i=0;i<=n;i++){
      fscanf(fi,"%lf%lf",&t1,&t2);
      *(x+i)=t1,*(y+i)=t2;
      if(t<fabs(t2))t=fabs(t2);
   }
/*计算三次样条插值函数所需参数*/
   for(i=1;i<=n;i++)
      *(h+i)=*(x+i)-*(x+i-1);
   for(i=1;i<=n-1;i++){
      *(u+i)=*(h+i)/(*(h+i)+*(h+i+1));
      *(v+i)=1-*(u+i);
      *(g+i)=6*((*(y+i+1)-*(y+i))/(*(h+i+1))-
           (*(y+i)-*(y+i-1))/(*(h+i)))/(*(h+i)+*(h+i+1));
   }
   RunAfterMethod(u,v,g,M);
   initgraph(&gd,&gm,"d:\\tc20\\bgi");
   line(0,240,640,240);
/*按比率缩放屏幕显示*/
   t1=600/(*(x+n)-*x),t2=100/t;
   for(i=0;i<=n;i++)
      circle(20+t1*(*(x+i)-*x),240-t2**(y+i),2);
   moveto(20,240-t2**y);
   for(t=*x;t<=*(x+n);t+=1/t1)
      lineto(20+t1*(t-*x),240-t2*S(t,x,y,h,M));
   free(x),free(y),free(h);
   free(u),free(v),free(g);
   fclose(fi);getch();closegraph();
}