na1006ac.cpp

数值分析的c语言几种算法实现

#include <stdio.h>

 

#define MAX_N 20

 

void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, 

double a[], double b[], double c[], double d[]);

 

double S( double t, double Fmax, 

            int n, double x[], double a[], double b[], double c[], double d[] );

 

int main()

{

    int n, Type, m, i;

    double x[MAX_N], f[MAX_N], a[MAX_N], b[MAX_N], c[MAX_N], d[MAX_N];

    double s0, sn, Fmax, t0, tm, h, t;

 

    while (scanf("%d", &n) != EOF) {

       for (i=0; i<=n; i++) 

           scanf("%lf", &x[i]);

       for (i=0; i<=n; i++) 

           scanf("%lf", &f[i]);

       scanf("%d %lf %lf %lf", &Type, &s0, &sn, &Fmax);

 

       Cubic_Spline(n, x, f, Type, s0, sn, a, b, c, d);

       for (i=0; i<=n; i++)

           printf("%12.8e %12.8e %12.8e %12.8e \n", a[i], b[i], c[i], d[i]);

 

       scanf("%lf %lf %d", &t0, &tm, &m);

       h = (tm-t0)/(double)m;

       for (i=0; i<=m; i++) {

           t = t0+h*(double)i;

           printf("f(%12.8e) = %12.8e\n", t, S(t, Fmax, n, x, a, b, c, d));

       }

       printf("\n");

    }

 

    return 0;

}


void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, 

double a[], double b[], double c[], double d[])
{
	double h[MAX_N];
	double aa[MAX_N];
	double ll[MAX_N];
	double u[MAX_N];
	double z[MAX_N];

	int i,j,k;
	for(i=0;i<=n;i++)
		a[i]=f[i];	
	for(i=0;i<n;i++)
		h[i]=x[i+1]-x[i];
	if(Type==2)
	{
	
		for(i=1;i<n;i++)
			aa[i]=3.0*((a[i+1]-a[i])/h[i]-(a[i]-a[i-1])/h[i-1]);
		ll[0]=1;
		u[0]=0;
		z[0]=0;

		for(i=1;i<n;i++)
		{
			ll[i]=2*(x[i+1]-x[i-1])-h[i-1]*u[i-1];
			u[i]=h[i]/ll[i];
			z[i]=(aa[i]-h[i-1]*z[i-1])/ll[i];
		}

		ll[n]=1;
		z[n]=0.0;
		c[n]=0.0;

		for(j=n-1;j>=0;j--)
		{
			c[j]=z[j]-u[j]*c[j+1];
			b[j]=(a[j+1]-a[j])/h[j]-h[j]*(c[j+1]+2*c[j])/3;
			d[j]=(c[j+1]-c[j])/(3*h[j]);
		}
		c[n]=sn/2.0;
		c[0]=s0/2.0;
	}
	else if(Type=1)
	{

		aa[0]=3*(a[1]-a[0])/h[0]-3*s0;
		aa[n]=3*sn-3*(a[n]-a[n-1])/h[n-1];
		for(i=1;i<n;i++)
			aa[i]=3.0*((a[i+1]-a[i])/h[i]-(a[i]-a[i-1])/h[i-1]);

		ll[0]=2*h[0];
		u[0]=0.5;
		z[0]=aa[0]/ll[0];

		for(i=1;i<n;i++)
		{
			ll[i]=2*(x[i+1]-x[i-1])-h[i-1]*u[i-1];
			u[i]=h[i]/ll[i];
			z[i]=(aa[i]-h[i-1]*z[i-1])/ll[i];
		}

		ll[n]=h[n-1]*(2-u[n-1]);
		z[n]=(aa[n]-h[n-1]*z[n-1])/ll[n];
		c[n]=z[n];

		for(j=n-1;j>=0;j--)
		{
			c[j]=z[j]-u[j]*c[j+1];
			b[j]=(a[j+1]-a[j])/h[j]-h[j]*(c[j+1]+2*c[j])/3;
			d[j]=(c[j+1]-c[j])/(3*h[j]);
		}

	}
	for(i=n;i>0;i--)
	{
		a[i]=a[i-1];
		b[i]=b[i-1];
		c[i]=c[i-1];
		d[i]=d[i-1];
	}
}

double S( double t, double Fmax, 

            int n, double x[], double a[], double b[], double c[], double d[] )
{
	if(t<x[0]||t>x[n]) return Fmax;
	int i,j,k;
	double temp;
	double xx;
	for(i=0;i<n;i++)
		if(t>=x[i]&&t<=x[i+1]) 
			break;
	i++;
	xx=a[i];
	temp=t-x[i-1];
	xx+=(b[i]*temp);
	temp*=(t-x[i-1]);
	xx+=(c[i]*temp);
	temp*=(t-x[i-1]);
	xx+=(d[i]*temp);
	return xx;
}