1007ac.txt

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

#include <stdio.h>

#include <math.h>

#define MAX_m 200

#define MAX_n 5



double f1(double x)

{

    return sin(x);

}



double f2(double x)

{

    return exp(x)*exp(x);

}



int OPA( double (*f)(double t), double x[], double w[], int m,

double c[], double *eps );



void print_results( int n, double c[], double eps)

{

    int i;



    printf("%d\n", n);

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

       printf("%8.4e ", c[i]);

    printf("\n");

    printf("error = %6.2e\n", eps);

    printf("\n");

}



int main()

{

    int m, i, n;

    double x[MAX_m], w[MAX_m], c[MAX_n+1], eps;



    m = 90;

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

       x[i] = 3.1415926535897932 * (double)(i+1) / 180.0;

       w[i] = 1.0;

    }

    eps = 0.001;

    n = OPA(f1, x, w, m, c, &eps);

    print_results(n, c, eps);



    m = 200;

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

       x[i] = 0.01*(double)i;

       w[i] = 1.0;

    }

    eps = 0.001;

    n = OPA(f2, x, w, m, c, &eps);

    print_results(n, c, eps);

	scanf("%d",&n);

    return 0;

}
#include <iostream>
//#include <cmath>
using namespace std;
int OPA( double (*f)(double t), double x[], double w[], int m,

double c[], double *eps )
{
    double p[MAX_n+1][MAX_n+1];
    double ss[MAX_n+1];
    double fz,fm;
    double temp[MAX_n+1];
    double xx;
    double err;
    int i,j,k;
    double p1,p2,pk;
    double s,a[MAX_n+1],b[MAX_n+1],aa[MAX_n+1];
	for(i=0;i<MAX_n+1;i++) c[i]=0;
	int flag=0;
	for(j=0;j<m-1;j++) {
		if(f(x[j])-f(x[j-1])!=0) 
		{
			flag=1;
			break;
		}
	}
	if( flag==0){ c[0]=f(x[0]);*eps=0.0;return 0;}
    for(i=0;i<MAX_n+1;i++)
        ss[i]=0;
    for(i=0;i<MAX_n+1;i++)
        for(j=0;j<MAX_n+1;j++)
            p[i][j]=0;
    for(i=0;i<MAX_n+1;i++){
        a[i]=b[i]=aa[i]=0;
    }
    p[0][0]=1;
    fz=0;
    fm=0;
    for(i=0;i<m;i++){
        fz+=(w[i]*f(x[i]));
        fm+=(w[i]);
    }
    aa[0]=fz/fm;
    fz=0;
    for(i=0;i<m;i++){
        fz+=(w[i]*x[i]);
    }
    a[1]=fz/fm;
    p[1][0]=-a[1];
    p[1][1]=1;
    fz=0;
    fm=0;
    for(i=0;i<m;i++){
        fz+=(w[i]*f(x[i])*(x[i]-a[1]));
        fm+=(w[i]*(x[i]-a[1])*(x[i]-a[1]));
    }
    aa[1]=fz/fm;
    ss[0]=aa[0]-aa[1]*a[1];
    ss[1]=aa[1];
    //for(i=0;i<6;i++)    cout << ss[i] << " " ; cout << endl;
    //cout << "c xishu " << ss[0] << " " << ss[1] << endl;
    err=0;
             for(i=0;i<m;i++)
             {
                pk=0;
                xx=1.0;
                for(j=0;j<MAX_n+1;j++){
                    pk+=(ss[j]*xx);
                    xx*=x[i];
                }
                err+=((f(x[i])-pk)*(f(x[i])-pk)*w[i]);
             }
            // cout << "error " << err << endl;
             if(err<=*eps){
                *eps=err;
                for(i=0;i<MAX_n+1;i++)
                 c[i]=ss[i];
                 return 1;
             }
    //cout << *eps << endl;
    for(k=2;k<MAX_n+1;k++)
    {
            fz=0;
            fm=0;
            for(i=0;i<m;i++){
                xx=1.0;
                p1=0;
                for(j=0;j<MAX_n+1;j++)
                {
                     p1+=(p[k-1][j]*xx);
                     xx*=x[i];
                }
                fz+=(w[i]*x[i]*p1*p1);
                fm+=(w[i]*p1*p1);
            }
            a[k]=fz/fm;
            fz=0;
            fm=0;
            for(i=0;i<m;i++){
                xx=1.0;
                p1=0;
                p2=0;
                for(j=0;j<MAX_n+1;j++)
                {
                     p1+=(p[k-1][j]*xx);
                     p2+=(p[k-2][j]*xx);
                     xx*=x[i];
                }
                fz+=(w[i]*p1*p1);
                fm+=(w[i]*p2*p2);
            }
            b[k-1]=fz/fm;
            for(i=0;i<MAX_n+1;i++)
             temp[i]=0;
            for(i=1;i<MAX_n+1;i++){
              temp[i]=p[k-1][i-1];
            }
             for(i=0;i<MAX_n+1;i++){
              temp[i]-=(p[k-1][i]*a[k]);
              temp[i]-=(p[k-2][i]*b[k-1]);
             }
             for(i=0;i<MAX_n+1;i++)
               p[k][i]=temp[i];

             fz=0.0;
             fm=0.0;
             for(i=0;i<m;i++){
                pk=0.0;
                xx=1.0;
                for(j=0;j<MAX_n+1;j++){
                    pk+=(p[k][j]*xx);
                    xx*=x[i];
                }
                fz+=(w[i]*f(x[i])*pk);
                fm+=(w[i]*pk*pk);
             }
             aa[k]=fz/fm;
             
             for(j=0;j<MAX_n+1;j++)
                ss[j]=ss[j]+p[k][j]*aa[k];
                
             err=0;
             for(i=0;i<m;i++)
             {
                pk=0;
                xx=1.0;
                for(j=0;j<MAX_n+1;j++){
                    pk+=(ss[j]*xx);
                    xx*=x[i];
                }
                err+=(f(x[i])-pk)*(f(x[i])-pk)*w[i];
             }
             if(k>=MAX_n||err<=*eps){
                *eps=err;
                for(i=0;i<MAX_n+1;i++)
                 c[i]=ss[i];
                 return k;
             }
       // for(i=0;i<6;i++)    cout << ss[i] << " " ; cout << endl;
    }
    /*for(i=0;i<6;i++)
        c[i]=ss[i];
    *eps=err;
    return 5;*/
}