scyt.cpp

求解差值真的范围场馆哟个自行车不是你沙发是但书

#include <iostream.h>
#include <fstream.h>

void input_scyt(int &n, double x[], double y[], double &yp1, double &ypn, double &x0)
{
	int i;
	ifstream infile("input_scyt.txt");
	infile>>n;
	for(i=1;i<=n;i++)
	{
		infile>>x[i];
	}
	for(i=1;i<=n;i++)
	{
		infile>>y[i];
	}
	infile>>yp1;
	infile>>ypn;
	infile>>x0;
}

double splint(double x[], double y[],  int n, double yp1, double ypn, double &x2)
{
	double u[100],y2[100],aaa,sig,p,bbb,ccc,qn,un,xa[100],ya[100],y2a[100];
	int i,k;
	for(i=1;i<=n;i++)
	{
		xa[i]=x[i];
		ya[i]=y[i];
	}
    if( yp1 > 9.9e+29 )
	{
        y2[1] = 0;
        u[1] = 0;
	}
    else
	{
        y2[1] = -0.5;
        aaa = (y[2] - y[1]) / (x[2] - x[1]);
        u[1] = (3.0 / (x[2] - x[1])) * (aaa - yp1);
    }
    for (i = 2; i<=n-1; i++)
	{
        sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
        p = sig * y2[i - 1] + 2.0;
        y2[i] = (sig - 1.0) / p;
        aaa = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
        bbb = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
        ccc = x[i + 1] - x[i - 1];
        u[i] = (6.0 * (aaa - bbb) / ccc - sig * u[i - 1]) / p;
    }
    if (ypn > 9.9e+29 )
	{
        qn = 0.0;
        un = 0.0;
	}
    else
	{
        qn = 0.5;
        aaa = ypn - (y[n] - y[n - 1]) / (x[n] - x[n - 1]);
        un = (3.0/ (x[n] - x[n - 1])) * aaa;
    }
    y2[n] = (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1.0);
    for (k = n - 1;k>=1;k--)
	{
        y2[k] = y2[k] * y2[k + 1] + u[k];
	}
	
	int klo,khi,k2;
	double h,a,b,aaa2,bbb2,yin;
    klo = 1;
    khi = n;
	for(i=1;i<=n;i++)
	{
		y2a[i]=y2[i];
	}
loop:   if (khi - klo > 1 )
	{
		k2 = (khi + klo) / 2;
		if (xa[k2] > x2)
			khi = k2;
		else
		{
			klo = k2;
		}
           goto loop;
	}
    h = xa[khi] - xa[klo];
    if (h == 0 )
	{
		cout<<"  pause  'bad  xa  input'"<<endl;
		return(0);
	}
	a = (xa[khi] - x2) / h;
	b = (x2 - xa[klo]) / h;
	aaa2 = a * ya[klo] + b * ya[khi];
	bbb2 = (a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi];
	yin = aaa2 + bbb2 * h*h /6.0;
	return(yin);
}

void main(void)
{
	double x[6],y[6],yp1,ypn,xin,yin;
	int n;
	input_scyt(n,x,y,yp1,ypn,xin);
	yin=splint(x,y,n,yp1,ypn,xin);
	cout<<yin<<endl;
}