/***********************************************************************
* 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 <stdio.h>
#include <math.h>

      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