/*********************************************************************** * This program solves a nonlinear system of two variables: * * * * f(x,y) = 0 * * g(x,y) = 0 * * * * using the Non_linear_system procedure. * * -------------------------------------------------------------------- * * REFERENCE: "Mathematiques en Turbo-Pascal part 1, By M. Ducamp and * * A. Reverchon, Editions EYROLLES, Paris, 1987" [BIBLI 03] * * * * C++ version by J-P Moreau, Paris * * (www.jpmoreau.fr) * * -------------------------------------------------------------------- * * SAMPLE RUN: * * (Solve the non linear system * * x^2 + x + y^2 - 2 = 0 * * x^2 + y - y^2 - 1 + ln(x) = 0 ) * * * * Desired precision........... : 1e-6 * * Maximal number of iterations : 100 * * Componants of starting vector: 1 1 * * * * Results in file nlinsyst.lst: * * (Error status: 0) * * * * -------------------------------------------------------------------- * * Method for nonlinear systems of two equations * * -------------------------------------------------------------------- * * * * System to be solved: * * x^2 + x + y^2 - 2 = 0 * * x^2 + y - y^2 - 1 + ln(x) = 0 * * * * Starting vector: * * 1.000000 1.000000 * * * * Error bound = 1e-6 * * Maximal number of iterations = 100 * * * * * * Solution vector: * * 0.915554 0.496191 * * * * Number of iterations: 6 * * -------------------------------------------------------------------- * * * ***********************************************************************/ #include #include int n, // size of system maxit, // maximal number of iterations it; // number of iterations performed double eps, // desired accuracy x0,yy0, // starting point x,y; // approximate solution FILE *fp; // output file int error; // error code: 0 = OK, 1 = error in evluating // f(x,y) or g(x,y), 2 = singular system double f(double x,double y) { return (x*x+x+y*y-2.0); } double g(double x,double y, int *rc) { *rc=0; if (x>0.0) return (x*x+y-y*y-1.0+log(x)); else { *rc=1; return 0; } } void Non_linear_system(double x0,double yy0, double prec, int maxiter, double *x, double *y, int *iter) { double h = 0.01; int rc; double a,b,c,d,t,m,n,p,q; error=1; *iter=0; *x=x0; *y=yy0; do { *iter=*iter+1; if (*iter > maxiter) return; a=f(*x+h,*y); b=g(*x+h,*y,&rc); if (rc!=0) return; a=(a-f(*x-h,*y))/2.0/h; b=(b-g(*x-h,*y,&rc))/2.0/h; if (rc!=0) return; c=f(*x,*y+h); d=g(*x,*y+h,&rc); if (rc!=0) return; c=(c-f(*x,*y-h))/2.0/h; d=(d-g(*x,*y-h,&rc))/2.0/h; if (rc!=0) return; t=a*d-b*c; if (fabs(t)<1e-12) { error=2; return; } m=f(*x,*y); n=g(*x,*y,&rc); if (rc!=0) return; p=(m*d-n*c)/t; q=(n*a-m*b)/t; *x=*x-p; *y=*y-q; } while ((fabs(p)+fabs(q)) > prec); error=0; } void main() { // -------------------- read input -------------------------------- n=2; // size of system printf("\n Desired precision..: "); scanf("%lf",&eps); printf(" Maximal number of iterations : "); scanf("%d",&maxit); printf(" Componants of starting vector:\n"); printf(" "); scanf(" %lf %lf",&x0,&yy0); // ------------ print input for checking purposes ----------------- fp = fopen("nlinsyst.lst","w"); fprintf(fp,"----------------------------------------------\n"); fprintf(fp," Method for nonlinear systems of two equations\n"); fprintf(fp,"----------------------------------------------\n\n"); fprintf(fp," System to be solved:\n"); fprintf(fp," x^2 + x + y^2 - 2 = 0\n"); fprintf(fp," x^2 + y - y^2 - 1 + ln(x) = 0\n"); fprintf(fp,"\n Starting vector:\n"); fprintf(fp," %f %f\n",x0,yy0); fprintf(fp,"\n Error bound = %f\n",eps); fprintf(fp," Maximal number of iterations = %d\n",maxit); // ------------ solve nonlinear system ---------------------------- Non_linear_system(x0,yy0,eps,maxit,&x,&y,&it); // --------------------- print solution --------------------------- fprintf(fp,"\n\n Solution vector:\n"); fprintf(fp," %f %f\n",x,y); fprintf(fp,"\n Number of iterations: %d\n", it); fprintf(fp,"----------------------------------------------\n\n"); fclose(fp); printf("\n Results in file nlinsyst.lst.\n"); printf(" Error status: %d\n\n",error); } // End of file nlinsyst.cpp