Re: Colnew Help

/* This is Example 4 with the notations as in  http://arxiv.org/pdf/hep-th/0407005 */

y:1.9d0$ aleft : 0.0d0$
aright : y$ ntol:4$
depends([p,f2,f3],t)$
lvfs:[t,p,f2,f3]$
list_of_eqns: [diff(p,t,2)= -0.5*(1/t+1/(t-1.0)+1/(t-y))*diff(p,t,1)+
((v-f2-f3)/t+f2/(t-1.0)+f3/(t-y))*p , diff(f2,t,1)=0, diff(f3,t,1)=0]$
list_of_conds:[(diff(p,t,1)-2.0*p*(v-f2-f3))(0)=0,
(p -1.0)(1)=0, (diff(p,t,1)-2.0*p*f2)(1)=0, (diff(p,t,1)-2.0*p*f3)(y)=0]$
list_of_guess:[2.0*t-1.0,1.0,1/(2.0*y-1.0)]$


load("colnew_help3.mac")$


use_result():= (xx:[0.2],zval: colnew_appsln(xx, mstar, fspace, ispace),
  print(zval))$ 
list_of_params: makelist(0.00001d0+0.2*k,k,0,9)$
for v in list_of_params do run_colnew(v)$

load(colnew)$
load("colnew_discov.lisp")$

colnew_symbs: [aleft,aright,lvfs,list_of_eqns,ntol,list_of_conds,list_of_guess]$

read_data(symb) ::= buildq([symb:symb], if (symb = 'symb) then
  (print("Define ", symb, ".Terminate by ;"), symb: read(" ",symb,": "))
  else symb)$

colnew_symbs: [aleft,aright,lvfs,list_of_eqns,ntol,list_of_conds,list_of_guess]$
for symb in colnew_symbs do (
if     (symb = aleft) then aleft:float(read_data(aleft))
elseif (symb = aright) then aright:float(read_data(aright))
elseif (symb = lvfs) then (
  print("Give the list of variable and functions in equations, like [x,F,G]"),
    lvfs:read_data(lvfs),
    depends(rest(lvfs),first(lvfs)))
elseif (symb = list_of_eqns) then (
  print("Define the list of the differential equations."),
  list_of_eqns: read_data(list_of_eqns),
  eqns_to_funs(list_of_eqns),
  define(funmake(df,cons(var0,var_list)), jacobian(apply(f,cons(var0,var_list)),var_list)))
elseif (symb = ntol) then (
  ntol:read_data(ntol),
  if (ntol < mstar) then print ("You need to edit ltol and tol")) 
elseif (symb = list_of_conds) then (
  print("Define the ", mstar, "boundary conditions."),
  list_of_conds: read_data(list_of_conds),
  conds_to_zconds(list_of_conds, aleft, aright, max_subst_list, m_list),
  fixpnt: float(fixpnt),
  ZETA:float(ZETA),
  define(funmake(dg, var_list), jacobian(apply(g,var_list),var_list)))
elseif (symb = list_of_guess) then (
  list_of_guess: read_data(list_of_guess),
  define(guess(var0), build_guess(list_of_guess,m_list,ncomp,var0))))$

build_ipar(m_list, mstar, ncomp, nfixpnt,ntol)$
fspace:makelist(0d0,k,1,ipar[5])$
ispace:makelist(0,k,1,ipar[6])$

print("ncomp = ", ncomp, ", m_list = ", m_list, ", mstar = ", mstar,
", ntol = ", ntol, ", max_subst_list = ", max_subst_list, "");
print("ZETA: ", ZETA);
print("ipar: ", ipar);
print("ltol: ", ltol);
print("tol: ",tol);
print("fixpnt: ",fixpnt);


iflag:1$
print("Please give floating numerical values to parameters.")$
ipar[7]:1$   /* to suppress abundant colnew printout */
break();
 
  
/* Assume param is free in f,g, etc but fixed in fsub,dfsub,etc. */

run_colnew(param):= (
define(funmake (fsub, cons(var0,var_list)), apply (f, cons(var0,var_list))),
define(funmake (dfsub, cons(var0,var_list)), apply(df,cons(var0,var_list))), 
define(funmake (gsub, cons(i,var_list)), buildq([gval:apply (g, var_list)], gval[i])),
define(funmake (dgsub, cons(i,var_list)), buildq([dgval:apply(dg,var_list)], row(dgval,i)[1])),
[iflag, fspace, ispace] : colnew_expert(ncomp, m_list, aleft, aright, ZETA,
  ipar, ltol, tol, fixpnt, ispace, fspace, iflag, fsub, dfsub, gsub,
  dgsub, guess),
use_result(),  /* must be defined in the main program */
ipar[9]:3,     /* run continuation after first step */
ipar[3]: ispace[1])$

;;; Input in the form:  diff(F,x,n)=... (lower order terms), conditions in the form
;;; C(F,...)(point)=0                     . Example (see example3):
;;; depends([G,H],x)$ 
;;; list-of-eqns:[ diff(G,x,2)=L^2*s*(G-1)-L*((3-n)/2*H*diff(G,x)+(n-1)*diff(H,x)*G),
;;; diff(H,x,3)= L^3*(1-G^2)+L^2*s*diff(H,x)-L*((3-n)/2*H*diff(H,x,2)+n*(diff(H,x))^2)]$
;;; list-of-conds:[G(0),H(0),diff(H,x)(0),(G-1)(1),diff(H,x)(1)]$
;;; Further implicit parameters, like L, s, n above may appear and used in continuation
;;; Output: n, variables z1,...,zn, and function f(x,z1,....,zn), a list of all eqs.
;;; Idem for conditions a list g(z1,....,zn) of all conds on cond points.

;;; The functions f(x,z1,....,zn) and g(z1,....,zn) are stored as maxima functions, hence
;;; define special variables $f, $g to which they are attached.
;;; Also give access to the many parameters of colnew: ncomp, mstar, ntol etc.

(declaim (special $f $g  $zvec $dmval $var0 $var_list $m_list $ntol $nfixpnt $fixpnt
		  $ncomp $mstar $aleft $aright $ipar |$zeta|))

(defun eq-disc (eqn)

  "Treat one maxima equation. It should be of the form diff(f,x,n)=... Check that."

  (if (not (equal (caar eqn) 'MEQUAL))
      (merror (intl:gettext "eqn:  needs an equation.")))
  (let* ((lhs (cadr eqn))
	 (diffop (caar lhs))
	 (functn (cadr lhs))
	 (variable (caddr lhs))
	 (order (cadddr lhs)))
    (if (not (equal diffop '%DERIVATIVE))
	(merror (intl:gettext "eqn: needs a differential equation.")))
    (if (> order 4) (merror (intl:gettext "Maximum degree of equation 4")))
    (values order functn variable)))

(defun cond-disc (condt)
  
  "Treat one cond, extract the condition and the point at which it is given. To keep things simple
   assume the condition is given in the form C(F,G,..)(p)=0 at point p. For example:
   :lisp #$$ (G+diff(G,x)-1)(3)=0$ -> ((MEQUAL) ((MQAPPLY) ((MPLUS) $g (($DIFF) $g $X) ((MMINUS) 1)) 3) 0)"
  
  (if (not (and (eql (caar condt) 'MEQUAL) (eql (caddr condt) 0)))
      (merror (intl:gettext "condt: must be of the form C(p)=0.")))
  (let* ((lhs (cadr condt)) 
	 (condit (if (eql (caar lhs) 'MQAPPLY) (cadr lhs) (caar lhs)))
	 (point (if (eql (caar lhs) 'MQAPPLY) (caddr lhs) (cadr lhs))))
    (values condit point)))


    
(defun $eqns_to_funs (list-of-eqns)

  "Treat a maxima list of maxima equations and a list of constraints. Outputs ncomp, mstar,
    m-list, for colnew, and assigns maxima function f and several maxima variables."

  (if (not (equal (caar list-of-eqns) 'MLIST))
      (merror (intl:gettext "Needs a list of equation.")))
  (let ((mstar 0) (ncomp 0) (index 0) m-list  newvar subst-list
	var0 var-list funs-list)
    (setq list-of-eqns (cdr list-of-eqns)) ; from maxima list to lisp list

;;; Treat equations
    (dolist (eqn list-of-eqns)
      (incf ncomp)
      (multiple-value-bind (order functn variable) (eq-disc eqn)
	(if (= 0 index) (setf var0 variable)
	    (if (not (equal variable var0))
		(merror (intl:gettext "Only one variable!"))))
	(push order m-list)
;;;  Introduce the Z-vector
	(dotimes (i order)
	  (setq newvar (intern (concatenate 'string "$Z" (write-to-string (incf index)))))
	  (push newvar var-list)
	  (if (= i 0)
	      (push `((MEQUAL)  ,functn ,newvar) subst-list)
	      (push `((MEQUAL)
		      ((%DERIVATIVE) ,functn ,variable ,i) ,newvar)
		    subst-list)))))
    (setf var-list (nreverse var-list))
    (if (> ncomp 20) (merror (intl:gettext "Maximum number of equations 20")))
    (setf mstar (reduce #'+ m-list))
    (setq m-list (cons '(MLIST) (nreverse m-list)))
    (if (> mstar 40) (merror (intl:gettext "Maximum sum of degrees 40")))
    (let ((max-subst-list (cons '(MLIST) subst-list)))
      (dolist (eqn list-of-eqns)
	; Doesn't work if max-subst-list is reversed!
	(push (simplify ($substitute max-subst-list (caddr eqn))) funs-list))
      (setq funs-list (nreverse funs-list))    
;;; Export variables to maxima
      (meval  `((MSETQ) $max_subst_list ,max-subst-list)))
    (meval  `((MSETQ) $mstar ,mstar))
    (meval  `((MSETQ) $ncomp ,ncomp))
    (meval  `((MSETQ) $m_list ,m-list))
    (meval  `((MSETQ) $var_list ,(cons '(MLIST) var-list)))
    (meval  `((MSETQ) $var0 ,var0))
    (meval  `((MDEFINE) ,(cons '($f) (list* var0  var-list)) ,(cons '(MLIST) funs-list))))
  'DONE)

(defun $conds_to_zconds(list-of-conds $aleft $aright $max_subst_list $m_list)
  
  "Express conditions in terms of the zvector, put them in order and discover if there are interior points to be
   inserted into the list of fix points."
  
  (setq list-of-conds (cdr list-of-conds))
  (if (not (eql (length list-of-conds) $mstar))
      (merror (intl:gettext "Need to have as many conditions as free parameters")))
  (let ((conds-list ()) (fixpnt ()) (zet ()) nfixpnt)
    (dolist (condt list-of-conds)
      (multiple-value-bind (condit point) (cond-disc condt)
	(push (list point (simplify ($substitute $max_subst_list  condit))) conds-list)))
    (setq conds-list (nreverse conds-list))
;;; Ensure points of ZETA are in order and constraints are sorted accordingly
    (setq conds-list (stable-sort conds-list #'< :key #'car))
    (setq zet (mapcar #'car conds-list))
    (setq conds-list (mapcar #'cadr conds-list))
    (dolist (zz zet) (if (and (> zz $aleft) (< zz $aright)
			      (not (member zz fixpnt)))
			 (push zz fixpnt)))
    (setq fixpnt (nreverse fixpnt))
    (setq nfixpnt (length fixpnt))
;;; Define  maxima function g, and maxima variables fixpnt ZETA (zeta is the function zeta)
    (meval  `((MDEFINE) ,(cons '($g) (cdr $var_list)) ,(cons '(MLIST) conds-list)))
;;; |$zeta| becomes ZETA in maxima (Manual, section 37.1) Other maxima variables
    (meval  `((MSETQ) |$zeta| ,(cons '(MLIST) zet)))
    (meval  `((MSETQ) $fixpnt ,(cons '(MLIST) fixpnt)))
    (meval  `((MSETQ) $nfixpnt ,nfixpnt)))
  'DONE)
    


(defun $build_guess ($list_of_guess $m_list $ncomp $var0)

  "list-of-guess is the init-guess, a maxima list of ncomp functions, giving a guess for the
   underivated functions. Construct z-vector and dmval in terms of var0, the variable of diff. eqns."

  (if (not (and (eql (caar $list_of_guess) 'MLIST)
		(eql ($length $list_of_guess) $ncomp)))
      (merror (intl:gettext "Guess must be a maxima list of functions.
                             of length ncomp")))
  (let (z-vec  dmval (ml (cdr $m_list)))
    (dolist (gfun  (cdr $list_of_guess))
      (let ((gexp gfun))
;;; For each guess compute derivatives of guess as indicated by m-list
	(dotimes (k (car ml))
	  (push gexp z-vec)
	  (setf gexp (simplify (meval `(($diff) ,gexp ,$var0)))))
	(push gexp dmval)
	(setq ml (cdr ml))))
    (setq z-vec (nreverse z-vec))
    (setq dmval (nreverse dmval))
    (setq $zvec (cons '(MLIST) z-vec))
    (setq $dmval (cons '(MLIST) dmval))
    (list '(MLIST) $zvec $dmval)))


;;; Example:  list_of_guess: [cos(x) + x^2,x^5]$ m_list:[2,3]$ ncomp:2$ var0:x$
;;; result: [[cos(x)+x^2,2*x-sin(x),x^5,5*x^4,20*x^3],[2-cos(x),60*x^2]]


(defun $build_ipar ($m_list $mstar $ncomp $nfixpnt $ntol &optional ($nmax 40) ($nsubint 5))
  
  ;; by default nmax=40 ntol=mstar tolerances set to 1d-5 ncomp=5
  ;; ntol is the number of tolerances, ltol specifies to what Z components they apply, tol their values.
  ;; for Problem 3 ipar=[1,4,10,2,40000,2500,1,0,1,0,0]
  ;; nmax is the maximum number of subintervals, ispace and fspace are resp  nmax*nsizei  nmax*nsizef.
  ;; One can add further fixed points in the array fixpnt and increase accordingly nfixpnt for example to
  ;; avoid a singularity in the equations, since the equations are never applied on fixed points. See Problem 4.
  
  (let ((mmax (apply 'max (cdr $m_list))))
    (if (= $ncomp 0)  (setq $ncomp (max (1+ mmax) (- 5 mmax)))
	(progn (if (< $ncomp mmax) (setq $ncomp mmax))
	       (if (> $ncomp 7) (setq $ncomp 7))))
    (if (= $nsubint 0) (setq $nsubint 5))
    (if (< $ntol 1) (setq $ntol 1))
    (if (> $ntol $mstar) (setq $ntol $mstar))
    (prepare-tol $ntol $mstar) ; fill tol and ltol with basic defaults
    (let* ((kd (* $nsubint $ncomp))
	   (kdm (+ kd $mstar))
	   (nsizei (+ 3 kdm))
	   (nsizef (+ 4 (* 3 $mstar) (* kdm (+ 5 kd)) (* 4 $mstar $mstar)))
	   ;; space is cheap one sets  (2*mstar-nrec)-> 2*mstar and also finally
	   (ispace (max (* $nmax nsizei) 2500))
	   (fspace (max (* $nmax nsizef) 40000))
	   ;; Here one sets ipar[9]=1 assuming there is a guess, for continuation one has to set
	   ;; ipar[9]=3 and ipar[3]:ispace[1] see the example in Problem 3.
	   (ipar (list '(MLIST) 1 $ncomp $nsubint $ntol fspace ispace 0 0 1 0 $nfixpnt)))

      (setq $ipar ipar)))
  'DONE)		


(defun prepare-tol ($ntol $mstar)
  
  ;; ltol is filled with increasing indices k (1<=k<=mstar) for which the tolerance has to be checked.
  ;; tol gives the corresponding tolerances. For example in problem 3 one has mstar=5, but ntol=2, and one chooses
  ;; ltol=[1,3], so Z1 and Z3 are checked, and tol=[1d-7,1d-7]. To see the correspondance with original variables
  ;; print max_subst_list, here only G and H are checked. If the tolerance on Z_k is eps, that means that
  ;; |Z_k - Z_k^0| < eps * Z_k^0 + eps where Z_k^0 is the exact solution and of course this difference is estimated
  ;; by looking at the rate of convergence of the Z_k after successive subdivisions.
  
  (defparameter $ltol  (cons '(MLIST) (loop for j below $ntol collect (1+ j))))
  (defparameter $tol (cons '(MLIST) (loop for f below $ntol collect 1.0D-5)))
  (if (< $ntol $mstar) (format t "You should edit the arrays ltol and tol ~%")
      (format t "The array tol is filled by default, you can edit tol ~%"))
  'DONE)

Loading /home/michel/.maxima/maxima-init.lisp
Loading /home/michel/.maxima/maxima-init.mac
Maxima 5.48.1 https://maxima.sourceforge.io
using Lisp SBCL 2.5.8
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) y:1.9d0$

(%i2) aleft : 0.0d0$

(%i3) aright : y$

(%i4) ntol:4$

(%i5) depends([p,f2,f3],t)$

(%i6) lvfs:[t,p,f2,f3]$

(%i7) list_of_eqns: [diff(p,t,2)= -0.5*(1/t+1/(t-1.0)+1/(t-y))*diff(p,t,1)+
((v-f2-f3)/t+f2/(t-1.0)+f3/(t-y))*p , diff(f2,t,1)=0, diff(f3,t,1)=0]$

(%i8) list_of_conds:[(diff(p,t,1)-2.0*p*(v-f2-f3))(0)=0,
(p -1.0)(1)=0, (diff(p,t,1)-2.0*p*f2)(1)=0, (diff(p,t,1)-2.0*p*f3)(y)=0]$

(%i9) list_of_guess:[2.0*t-1.0,1.0,1/(2.0*y-1.0)]$

(%i10) load("colnew_help3.mac")$
Give the list of variable and functions in equations, like [x,F,G] 
Define the list of the differential equations. 
Define the  4 boundary conditions. 
The array tol is filled by default, you can edit tol 
ncomp =  3 , m_list =  [2, 1, 1] , mstar =  4 , ntol =  4 , max_subst_list =  
                                                             dp
                                          [f3 = z4, f2 = z3, ── = z2, p = z1]  
                                                             dt
ZETA:  [0.0, 1.0, 1.0, 1.9] 
ipar:  [1, 3, 5, 4, 40000, 2500, 0, 0, 1, 0, 1] 
ltol:  [1, 2, 3, 4] 
tol:  [1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5] 
fixpnt:  [1.0] 
Please give floating numerical values to parameters. 

Entering a Maxima break point. Type 'exit;' to resume.
_exit;
(%i11) use_result():= (xx:[0.2],zval: colnew_appsln(xx, mstar, fspace, ispace),print(zval))$

(%i12) list_of_params: makelist(0.00001d0+0.2*k,k,0,9)$

(%i13) for v in list_of_params do run_colnew(v)$
[[- 0.37440899999619803, 1.7180073231152533, 0.8590075884393673, 
                                                          0.3373678645255313]] 
[[- 0.4181243778331842, 1.692717871309146, 0.9270126805146587, 
                                                          0.3714363278793304]] 
[[- 0.46516538992082285, 1.668234722091104, 1.0002886037429792, 
                                                          0.4051263445182374]] 
[[- 0.515142209185673, 1.643281449790386, 1.0791425591666848, 
                                                         0.43844377727383627]] 
[[- 0.5674149323450435, 1.6163968495521828, 1.163697787796863, 
                                                          0.4713601621913387]] 
[[- 0.6211324864560417, 1.5860537083169028, 1.2538759260131196, 
                                                          0.5038210392323347]] 
[[- 0.6753085929900787, 1.5508002653627255, 1.3494108272733665, 
                                                          0.5357580123658756]] 
[[- 0.7289152458649232, 1.5093908954492372, 1.4498880436819064, 
                                                          0.5671012806913346]] 
[[- 0.7809707724456075, 1.460878882531993, 1.554796899330325, 
                                                          0.5977896540847439]] 
[[- 0.8306060253095319, 1.404660134412408, 1.663581618812131, 
                                                          0.6277765342369905]] 
(%i14)