d9r10.cpp

来自「工程算法 这是一个很有用的工程数值算法集锦」· C++ 代码 · 共 67 行

#include "iostream.h"
#include "math.h"
#include "stdlib.h"

void main()
{
    //program d9r10
    //driver for routine medfit
    int i,npt = 100;
    double a,b,siga,sigb,chi2,q,abdev,spread = 0.1;
    double x[101], y[101], sig[101];
    long idum = -1984;
    for (i = 1; i<=npt; i++)
	{
        x[i] = 0.1 * i;
        y[i] = -2.0 * x[i] + 1.0 + spread * gasdev(idum);
        sig[i] = spread;
    }
    int mwt = 1;
    fit(x, y, npt, sig, mwt, a, b, siga, sigb, chi2, q);
	cout.setf(ios::fixed|ios::left);
	cout.precision(4);
    cout<<"According to routine fit the result is:"<<endl;
    cout<<endl;
    cout<<"a = ";
	cout.width(12);
	cout<<a;
	cout<<"Uncertainty: ";
	cout.width(12);
    cout<<siga<<endl;
	cout<<endl;
    cout<<"b = ";
	cout.width(12);
	cout<<b;
	cout<<"Uncertainty: ";
    cout<<sigb<<endl;
    cout<<endl;
    cout<<"Chi-squared: ";
	cout.width(12);
	cout<<chi2<<" for   ";
	cout.width(5);
	cout<<npt<<"points";
    cout<<endl;cout<<endl;
    cout<<"Goodness-of-fit: ";
	cout.width(12);
	cout<<q;
    cout<<endl;cout<<endl;
    cout<<"According to routine medfit the result is:";
    cout<<endl;cout<<endl;
    medfit(x, y, npt, a, b, abdev);
    cout<<"a = ";
	cout.width(12);
	cout<<a;
    cout<<"b = ";
	cout.width(12);
	cout<<b;
    cout<<endl;cout<<endl;
    cout<<"Absolute deviation (per data point): ";
	cout.width(12);
	cout<<abdev;
    cout<<endl;cout<<endl;
    cout<<"(note: Gaussian spread is  ";
	cout<<spread;
	cout<<")";
	cout<<endl;cout<<endl;
}