'************************************************************************
'* 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 or 1991"    *
'*             [BIBLI 03].                                              *
'*                                                                      *
'*                                  Basic 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 + log(x) = 0 )                    *
'*                                                                      *
'* Desired precision........... : 1e-10                                 *
'* Maximal number of iterations : 30                                    *
'* 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 + log(x) = 0                                      *
'*                                                                      *
'*  Starting vector:                                                    *
'*  1.000000  1.000000                                                  *
'*                                                                      *
'*  Error bound = 1e-10                                                 *
'*  Maximal number of iterations = 30                                   *
'*                                                                      *
'*                                                                      *
'*  Solution vector:                                                    *
'*  .9155544492454276    .4961910934523857                              *
'*                                                                      *
'*  Number of iterations: 7                                             *
'* -------------------------------------------------------------------- *
'*                                                                      *
'************************************************************************
defint i-n
defdbl a-h,o-z

'integer n              size of system
'        maxit          maximal number of iterations
'        iter           number of iterations performed

'double  eps            desired accuracy
'        x0,y0          starting point
'        x,y            approximate solution

'integer ierror         error code:  0 = OK
'                                    1 = error in evluating f(x,y) or g(x,y)
'                                    2 = singular system

' begin main program
' -------------------- read input --------------------------------
  n=2    'size of system
  cls
  print
  input " Desired precision..: ", eps
  input " Maximal number of iterations : ", maxit
  print " Componants of starting point (x0,y0):"
  input "    x0 = ", x0
  input "    y0 = ", y0

' ------------ print input for checking purposes -----------------
  open "nlinsyst.lst" for output as #1
  print #1,"----------------------------------------------"
  print #1," Method for nonlinear systems of two equations"
  print #1,"----------------------------------------------"
  print #1,
  print #1," System to be solved:"
  print #1," x^2 + x + y^2 - 2          = 0"
  print #1," x^2 + y - y^2 - 1 + log(x) = 0"
  print #1,
  print #1," Starting vector:"
  print #1,"  ";x0;"  ";y0
  print #1,
  print #1," Error bound = "; eps
  print #1," Maximal number of iterations = "; maxit

' ------------ solve nonlinear system ----------------------------
' call Non_linear_system subroutine
  gosub 2000
' --------------------- print solution ---------------------------
  print #1,
  print #1,
  print #1," Solution vector:"
  print #1,"  ";x;"  ";y
  print #1,
  print #1," Number of iterations: "; iter
  print #1,"----------------------------------------------"
  print #1,
  close #1

  print
  print " Results in file nlinsyst.lst."
  print " Error status: "; ierror
  print

end 'of main program


1000 'Function f(xx,yy)
  f=xx*xx+xx+yy*yy-2#
return

1200 'Function g(xx,yy)
  irc=0
  if xx>0# then
    g=xx*xx+yy-yy*yy-1#+LOG(xx)
  else
    irc=1 : g=0#
  end if
return


2000 'Subroutine Non_linear_system
'double a,b,c,d,t,xm,xn,p,q
  h = 0.01#
  ierror=1 : iter=0
  x=x0 : y=y0
2100 iter=iter+1
  if iter > maxit then return
  xx=x+h : yy=y : gosub 1000
  a=f
  gosub 1200 : if irc<>0 then return
  b=g
  xx=x-h : gosub 1000
  a=(a-f)/2#/h
  gosub 1200 : if irc<>0 then return
  b=(b-g)/2#/h
  xx=x : yy=y+h : gosub 1000
  c=f
  gosub 1200 : if irc<>0 then return
  d=g
  yy=y-h : gosub 1000
  c=(c-f)/2#/h
  gosub 1200 : if irc<>0 then return
  d=(d-g)/2#/h
  t=a*d-b*c
  if abs(t)<1e-12 then
    ierror=2
    return
  end if
  xx=x : yy=y : gosub 1000
  xm=f
  gosub 1200 : if irc<>0 then return
  xn=g
  p=(xm*d-xn*c)/t
  q=(xn*a-xm*b)/t
  x=x-p : y=y-q
  if abs(p)+abs(q) > eps then goto 2100
  ierror=0
return

'End of file nlinsyst.bas