share/integequations/inteqn.mac

   1 eval_when([translate,batch,demo,load,loadfile],
   2 dv(a)::=buildq([a],define_variable(a,'a,any)),
   3 dvl(a)::=buildq([a],define_variable(a,[],list)),ttyoff:true)$
   4
   5 /*  (c) Copyright 1981 Massachusetts Institute of Technology  */
   6
   7 /* dynamalloc:true$ */
   8 define_variable(ieqnprint,true,boolean)$
   9 define_variable(techlist,'[first,all,vlfrnk,transform,collocate,
  10                 flfrnk2nd,tailor,fredseries,neumann,
  11                 flfrnk1st,abel,firstkindseries],list)$
  12 dv(tech)$
  13 dv(px)$
  14 dv(eqn)$
  15 dv(fx)$
  16 dv(uvar)$
  17 dv(xvar)$
  18 dv(kxu)$
  19 dv(ax)$
  20 dv(bx)$
  21 dv(iesoln)$
  22 dv(napprox)$
  23 dv(pxlist)$
  24 dv(eqnlist)$
  25 dv(inteqn)$
  26 dv(lowlim)$
  27 dv(pofx)$
  28 dv(iclist)$
  29 dvl(unklist)$
  30 dv(highlim)$
  31 dv(%c)$
  32 dv(eqno)$
  33 dvl(ueqnlist)$
  34 dv(x)$
  35 dv(xhat)$
  36 dv(opr)$
  37 dv(guesslist)$
  38 dv(guess)$
  39
  40 ieqn([arglist]):=
  41    if length(arglist)<2 then error("ieqn requires at least two arguments") else
  42         ieqn1(arglist[1], arglist[2],
  43                 if length(arglist)>2 then arglist[3] else [],
  44                 if length(arglist)>3 then arglist[4] else [],
  45                 if length(arglist)>4 then arglist[5] else [])$
  46
  47 opr(e):=if mapatom(e) then 'none else inpart(e,0)$
  48 arg(e):=inpart(e,1)$
  49 alias(exitfor, return)$
  50
  51 /* identifies functions unknown to macsyma */
  52 unfun(e):=is(?getchar(e,1)='?\$)$
  53
  54 /* ***** main program *****/
  55 ieqn1(eqnlist,pxlist,tech,napprox,guesslist):=
  56    block([iesoln, xvar, uvar, ax, bx, kxu, fx, eqn, px, guess, system,
  57           inflag, dispflag, singsolve, solveradcan, keepfloat],
  58         inflag:true, dispflag:false, system:false,
  59         if not listp(eqnlist) then eqnlist:[eqnlist],
  60         if not listp(pxlist) then pxlist:[pxlist],
  61         if length(eqnlist)#length(pxlist) then
  62           error("number of equations # number of unknowns"),
  63         for unk in pxlist do
  64           if mapatom(unk) or length(unk)#1 or not atom(xvar:arg(unk)) or
  65              numberp(xvar) then error(unk," improperly specified unknown"),
  66         if tech=[] then myprint("default 3rd arg, technique: ",tech:'first) else
  67         if not member(tech,techlist) then error(tech," invalid technique - see techlist"),
  68         if napprox=[] then myprint("default 4th arg, number of iterations or coll. parms.: ",napprox:1) else
  69         if not integerp(napprox) or napprox<=0 then error("napprox must be a pos. integer"),
  70         if guesslist=[] then myprint("default 5th arg, initial guess: ",guesslist:'none),
  71         if guesslist='none then guess:'none else
  72           (if not listp(guesslist) then guesslist:[guess:guesslist],
  73            if length(guesslist)#length(pxlist) then error("number of guesses # number of unknowns")),
  74         singsolve:solveradcan:keepfloat:true, iesoln:[],
  75         eqnlist:maplist(lambda([eqn],num(ratsimp(lhs(eqn)-rhs(eqn)))),eqnlist),
  76         if length(pxlist)=1 then (px:first(pxlist), eqn:first(eqnlist))
  77           else system:true,
  78         if not system then
  79            (try('vlfrnk), if tech='done then return(reverse(iesoln))),
  80         try('transform),
  81         if tech='done then return(reverse(iesoln)),
  82         for eqnlist in solvesys(eqnlist,pxlist) do
  83            (try('flfrnk2nd), try('tailor),
  84            if not system then try('fredseries), try('neumann)),
  85         if not system and firstkind() then
  86            (try('flfrnk1st), try('abel), try('firstkindseries)),
  87         try('collocate), return(reverse(iesoln)))$
  88
  89 /* invokes solution method catching and containing errors.  globals: tech */
  90 try(routine):=if member(tech,['all, 'first, routine]) then
  91         block([attempt], attempt:errcatch(catch(apply(routine,[]))),
  92           if attempt=[] then (?princ('? in\ ), ?princ(?stripdollar(routine)),
  93                 ?terpri()) else
  94             if first(attempt) and tech='first then tech:'done)$
  95
  96
  97 /* test for 1st kind and extracts parts.
  98   globals: eqn, kxu, uvar, xvar, ax, bx */
  99 firstkind():=block([integral, xpoints, consistent],
 100         if not freeof(px,eqn) then return(false),
 101         integral:catch(findint(eqn)),
 102         if integral=false then return(false),
 103         fx:mysolve(eqn,integral), if fx=[] then return(false),
 104         fx:first(fx), uvar:inpart(integral,2),
 105         kxu:lin(arg(integral), subst(uvar,xvar,px)),
 106         if kxu=false then return(false),
 107         ax:inpart(integral,3), bx:inpart(integral,4),
 108         if kxu[2]#0 then fx:ratsimp(fx-myint(kxu[2],uvar,ax,bx)),
 109         xpoints:mysolve(ax-bx,xvar), kxu:kxu[1],
 110         if xpoints=[] then return(true) else consistent:false,
 111         for xv in xpoints do
 112           if ratsubst(xv,xvar,fx)=0 then exitfor(consistent:true),
 113         return(consistent))$
 114
 115 /* returns first integral found in expr */
 116 findint(expr):=if not mapatom(expr) then
 117                 (if opr(expr)#'?%integrate then
 118                         (for i in expr do findint(i), false)
 119                  else throw(expr))$
 120
 121 /* globals: ieqnprint, iesoln */
 122 printsoln(soln,tech,order,kind):=block([dispflag, den],
 123         dispflag:ieqnprint, if kind#[] then kind:[kind],
 124         if order#[] then kind:cons(order,kind),
 125         if not listp(soln) then soln:[soln],
 126         soln:maplist(lambda([elmt], block([elem,den],
 127         elem:ratsimp(elmt), den:denom(elem),
 128         if numberp(den) and den#1 and opr(num(elem))="+" then
 129           multthru(elem) else elem)),soln),
 130         if rest(soln)=[] then soln:first(soln),
 131         iesoln:cons(?displine(cons(soln,cons(tech,kind))),iesoln),
 132         return(true))$
 133
 134 alias(myint, integrate)$
 135
 136 myprint(msg1,msg2):=if ieqnprint then print(msg1,msg2)$
 137
 138 /* if expr = a*var+b returns [a,b] else false */
 139 lin(expr,var):= block([lc],
 140         lc:ratsimp(diff(expr,var)),
 141         if freeof(var,lc) then [lc, ratsubst(0,var,expr)])$
 142
 143 /* solves eqn for unk and returns a list of actual solution values */
 144 mysolve(eqn,unk):=block([result], result:[],
 145         eqn:solve(eqn,unk), eqn:maplist('rhs, apply('ev,[eqn])),
 146         for ans in eqn do
 147           if freeof(unk, ans) then result:cons(ans, result),
 148         return(result))$
 149
 150 /* similar to mysolve but for systems. */
 151 solvesys(eqns,unks):=
 152         (eqns:solve(eqns,unks), eqns:apply('ev,[eqns]), maplist(lambda([el],
 153          if listp(el) then maplist('rhs,el) else [rhs(el)]), eqns))$
 154
 155 /* if expr can be expressed as sum(f[i](x)*g[i](u),i,1,n) then
 156   result is [[f1,g1],...[fn,gn]] else [].  globals: xvar  */
 157 sumfactors(expr,uvar):=block([other, rem],
 158         expr:catch(frnk(expr)), if expr=false then return([]),
 159         if freeof(uvar,expr) then return([[expr,1]]),
 160         if freeof(xvar,expr) then return([[1,expr]]),
 161         expr:expand(expr),
 162         other:mypartition(expr,uvar), expr:mypartition(expr,xvar),
 163         if other[1]<expr[1] then expr:other,
 164         rem:if expr[3]=0 then [] else [[expr[3],1]],
 165         if expr[4]#0 then rem:cons([1,expr[4]],rem),
 166         if expr[2]#0 then for fc in expr[2] do rem:cons(partition(fc,uvar),rem),
 167         return(rembox(rem)))$
 168
 169 /* converts e to finiterank form if possible otherwise throws false.
 170   globals: xvar, uvar  */
 171 frnk(e):=
 172         if freeof(xvar,e) or freeof(uvar,e) then e else
 173         block([op,ar,up], op:opr(e), ar:arg(e),
 174           if member(op,["+","*"]) then return(map('frnk,e)),
 175           if op="^" then (up:inpart(e,2),
 176            if integerp(up) then
 177                 if up>0 then return(frnk(ar)^up) else throw(false),
 178            up:plussplit(expand(up)),
 179            if up#[] and freeof(uvar,xvar,ar) then
 180                return(box(ar^up[1])*box(ar^up[2])) else throw(false)),
 181           if member(opr(ar),["*","+"]) then (e:partition(ar,xvar),
 182             if not freeof(uvar,e[2]) then throw(false)),
 183           if op='?%log and opr(ar)="*" then log(e[1])+log(e[2]) else
 184           if opr(ar)="+" then
 185             if op='?%sin  then sin(e[1])*cos(e[2])+cos(e[1])*sin(e[2]) else
 186             if op='?%cos  then cos(e[1])*cos(e[2])-sin(e[1])*sin(e[2]) else
 187             if op='?%sinh then sinh(e[1])*cosh(e[2])+cosh(e[1])*sinh(e[2]) else
 188             if op='?%cosh then cosh(e[1])*cosh(e[2])+sinh(e[1])*sinh(e[2]) else
 189             throw(false) else
 190           throw(false))$
 191
 192 /* if e=f(x)+g(u) returns [f(x),g(u)] else [].  globals: xvar, uvar */
 193 plussplit(e):=
 194         if opr(e)="+" then (e:partition(e,uvar),
 195             if freeof(xvar,e[2]) then e else []) else
 196         if freeof(uvar,e) then [e,0] else
 197         if freeof(xvar,e) then [0,e] else []$
 198
 199 /* solves exs for unks and substitutes in form */
 200 solveandsubst(exs,unks,form,tech):=block([soln,trial],
 201         if (soln:solve(exs,unks))=[] then
 202           (print("for a ",tech," solution substitute in the expression:"),
 203            ldisp(form), apply('print,cons("the values of ",unks)),
 204            print("that make the following expressions simultaneously zero"),
 205            apply('ldisp,exs), return(false)),
 206         for sol in apply('ev,[soln]) do (trial:apply('ev,[form,sol]),
 207           apply('printsoln,append([trial,tech],
 208             if tech='collocate then [napprox,testsoln(trial)] else [[],[]]))),
 209         return(true))$
 210
 211 /* tests solution for exactness.  globals: eqnlist, pxlist */
 212 testsoln(resultlist):=block([flag], apply('local,maplist('opr,pxlist)),
 213         flag:[], maplist('define,pxlist,resultlist),
 214         resultlist:apply('ev,[eqnlist,'integrate,'ratsimp]),
 215         for val in resultlist do
 216           if val#0 then exitfor(flag:'approximate),
 217         return(flag))$
 218
 219
 220 /* fout=h(x,u)+q(x)+r(u).  returns [n,h(x,u),q(x),r(u)] where
 221   n is the rank of fout.  globals: xvar, uvar */
 222 mypartition(fout,var):=block([qx, ru, rem, con],
 223         if opr(fout)#"+" or opr(fout:factorout(fout,var))#"+" then
 224                 return([1,[fout],0,0]),
 225         rem:con:qx:ru:0,
 226         for fc in fout do
 227           if freeof(uvar,fc) then
 228                 (if freeof(xvar,fc) then con:con+fc else qx:qx+fc) else
 229           if freeof(xvar,fc) then ru:ru+fc else rem:rem+fc,
 230         fout:if rem=0 then 0 else
 231              if opr(rem)="+" then length(rem) else (rem:[rem], 1),
 232         if qx#0 then (fout:fout+1, qx:qx+con, con:0),
 233         if ru#0 then (fout:fout+1, ru:ru+con),
 234         return([fout, rem, qx, ru]))$
 235
 236 /* replaces all integrals and unknown functions in ex with zero  */
 237 zeroint(ex):=
 238         if mapatom(ex) then ex else
 239         if opr(ex)='?%integrate or unfun(opr(ex)) then 0 else
 240         map('zeroint,ex)$
 241
 242 /* ----------- the following functions apply to 1st or 2nd kind eqns */
 243
 244 /* variable-limit finite-rank 1st or 2nd kind. reduction to ode.
 245   globals: eqn, xvar, px */
 246 vlfrnk():=block([initial, lowlim, unklist, inteqn, pofx, kind,
 247                         firstkind, difflist, iclist],
 248         pofx:px, initial:true, unklist:difflist:[], kind:'incomplete,
 249         firstkind:is(freeof(px,eqn)), inteqn:vconvert(eqn),
 250         iclist:[xvar=lowlim],
 251         for count thru length(unklist) do
 252           ( if firstkind then firstkind:false else
 253               (if initial then getics(), pofx:diff(pofx,xvar)),
 254             difflist:cons(inteqn,difflist), inteqn:diff(inteqn,xvar),
 255             if freeof('?%integrate, inteqn) then exitfor(false)),
 256         apply('remove,[opr(px), 'atvalue]),
 257         if not freeof('?%integrate, inteqn) then
 258           ( difflist:solve(difflist, unklist),
 259             if difflist=[] then throw(false),
 260  inteqn:apply('ev,[inteqn,difflist])),
 261         lowlim:derivdegree(inteqn,px,xvar), iclist:reverse(iclist),
 262         if lowlim=0 then (iclist:[], lowlim:3,
 263           if (pofx:mysolve(inteqn,px))#[] then (inteqn:first(pofx), kind:[])),
 264         if lowlim>2 or (pofx:ode2(inteqn,px,xvar))=false then
 265             return(printsoln(inteqn, 'vlfrnk, iclist, kind)),
 266         pofx:ratsimp(pofx),
 267         if initial and length(iclist)-1=lowlim and (inteqn:errcatch(
 268            apply(if lowlim=1 then 'ic1 else 'ic2, cons(pofx,iclist))))#[] then
 269              (pofx:first(inteqn), iclist:[]),
 270         if opr(pofx)="=" and lhs(pofx)=px then (pofx:rhs(pofx), kind:[]),
 271         if iclist=[] and kind#[] then iclist:0,
 272         printsoln(pofx, 'vlfrnk, iclist, kind))$
 273
 274 /* obtains initial conditions.
 275    globals: inteqn, xvar, lowlim, pofx, iclist, initial */
 276 getics():=block([val, init],
 277         val:at(inteqn, xvar=lowlim), init:at(pofx, xvar=lowlim),
 278         init:mysolve(val, init), if init#[] then
 279           ( init:first(init), atvalue(pofx, xvar=lowlim, init),
 280             iclist:cons(pofx=init, iclist)) else initial:false)$
 281
 282
 283 /* converts integrands in fun to finiterank form. upper limit must
 284 be xvar and lower limit must be a constant.
 285   globals: initial, lowlim, unklist, xvar */
 286 vconvert(fun):=
 287         if mapatom(fun) then fun else
 288         if opr(fun)#'?%integrate then map('vconvert, fun) else
 289         if not freeof(xvar,inpart(fun,3)) or
 290            inpart(fun,4)#xvar then throw(false) else
 291         block([intgr, newfun, int],
 292                 if lowlim='lowlim then lowlim:inpart(fun,3) else
 293                 if lowlim#inpart(fun,3) then initial:false,
 294                 intgr:sumfactors(arg(fun), inpart(fun,2)),
 295                 if intgr=[] then throw(false), newfun:0,
 296                 for term in intgr do
 297                   (int:substinpart(term[2], fun, 1),
 298                    if not member(int,unklist) then unklist:cons(int, unklist),
 299                    newfun:newfun+term[1]*int),
 300                 return(newfun))$
 301
 302
 303 /* laplace transform.   globals: eqnlist, pxlist, xvar  */
 304 transform():=block([teqnlist, translist, %s, ps, flag],
 305         ps:lambda([fun], laplace(fun, xvar, %s)), flag:false,
 306         teqnlist:maplist('ps,eqnlist), translist:maplist('ps,pxlist),
 307         for soln in solvesys(teqnlist,translist) do
 308          if freeof('?%integrate, '?%laplace, soln) then
 309            (soln:maplist(lambda([fun], ilt(fun,%s,xvar)), soln),
 310            ps:freeof('?%ilt, soln),
 311            printsoln(soln, 'transform, if ps then [] else 0, if ps then [] else 'incomplete),
 312            flag:flag or ps),
 313         return(flag))$
 314
 315 /* collocation.  globals: eqnlist, pxlist, napprox, xvar */
 316 collocate():=block([lowlim, highlim, unklist, %c, elist, form, incr, point,
 317                         name, listeqns, solnlist],
 318         apply('local, cons('%c,maplist('opr,pxlist))), lowlim:'minf,
 319 highlim:'inf,
 320         map('getlimits, eqnlist), if highlim<=lowlim then throw(false),
 321         solnlist:listeqns:unklist:[],
 322         for unk in pxlist do
 323         ( name:opr(unk), form:approx(name,napprox,xvar),
 324           apply('define,[unk,form]), solnlist:cons(form,solnlist),
 325           for parm:0 thru napprox-1 do unklist:cons(%c[name,parm],unklist)),
 326         elist:apply('ev,[eqnlist,'integrate,'expand]),
 327         if not freeof('?%integrate,elist) then throw(false),
 328         if freeof(xvar,elist) then listeqns:elist else (point:lowlim,
 329          incr:if napprox>1 then (highlim-lowlim)/(napprox-1) else 0,
 330          for i thru napprox do
 331            (listeqns:append(subst(point,xvar,elist),listeqns),
 332            point:point+incr)),
 333         solveandsubst(listeqns,unklist,reverse(solnlist),'collocate))$
 334
 335 /* obtains largest lower limit and least upper limit for collocation points
 336         globals: lowlim, highlim */
 337 getlimits(expr):=
 338         if not mapatom(expr) then
 339           if opr(expr)#'?%integrate then
 340                 for sub in expr do getlimits(sub) else
 341           block([low,high], low:inpart(expr,3), high:inpart(expr,4),
 342                 if not numberp(low) or not numberp(high) then throw(false),
 343                 lowlim:max(low,lowlim), highlim:min(high,highlim))$
 344
 345 /* approximation function for unknown function fun(var)  */
 346 approx(fun,nparms,var):=sum(%c[fun,i-1]*var^(i-1),i,1,nparms)$
 347
 348 /* ---- the following functions apply only to 2nd kind eqns */
 349
 350 /*  fixed-limit, finite-rank, 2nd kind  globals:  eqnlist, pxlist  */
 351 flfrnk2nd():=block([frlist, unklist, ueqnlist, eqno, %c],
 352         apply('local, cons('%c,maplist('opr, pxlist))),
 353         unklist:ueqnlist:[], eqno:0,
 354         frlist:maplist('fconvert, eqnlist),
 355         maplist('define, pxlist, frlist),
 356         ueqnlist:apply('ev,[ueqnlist,'integrate]),
 357         solveandsubst(ueqnlist, unklist, frlist, 'flfrnk2nd))$
 358
 359 /*  returns result of replacing all occurrences of
 360         'integrate(f(x,u),u,a,b) in fun with
 361         sum(q[j](x)*%c[j],j,1,n) where %c[j] is 'integrate(r[j](u),u,a,b)
 362         after expressing integrands in finite-rank form
 363         (here f(x,u) becomes sum(q[j](x)*r[j](u),j,1,n).
 364         globals: eqno (current number of subscripted unknowns %c),
 365                  ueqnlist (list of equations relating %c's)
 366                  unklist (list of all %c unknowns to be solved for)
 367                  xvar (independent variable) */
 368 fconvert(fun):=
 369         if mapatom(fun) then fun else
 370         if opr(fun)#'?%integrate then map('fconvert, fun) else
 371         if not freeof(xvar,inpart(fun,3)) or not freeof(xvar, inpart(fun,4))
 372          then throw(false) else
 373         block([newfun, intgrnd],
 374           intgrnd:sumfactors(arg(fun), inpart(fun,2)),
 375           if intgrnd=[] then throw(false),  newfun:0,
 376           for term in intgrnd do
 377                 (eqno:eqno+1, newfun:newfun+%c[eqno]*term[1],
 378                  ueqnlist:cons(%c[eqno]-substinpart(term[2],fun,1), ueqnlist),
 379                  unklist:cons(%c[eqno], unklist)),
 380           return(newfun))$
 381
 382 /*  taylor series. globals: eqnlist, pxlist, xvar, napprox  */
 383 tailor():=block([xhat, ufun, eqtn, tfun, neqns, vlist, order, value, fact],
 384         apply('local,append([ufun,eqtn,tfun],maplist('opr,pxlist))),
 385         map('getxhat, eqnlist), neqns:0, vlist:pxlist,
 386         for expr in eqnlist do
 387         ( neqns:neqns+1, eqtn[neqns]:expr, ufun[neqns]:first(vlist),
 388           vlist:rest(vlist), tfun[neqns]:value:subst(xhat,xvar,expr),
 389           atvalue(ufun[neqns], xvar=xhat, value)),
 390         fact:1, order:napprox,
 391         for i thru napprox do
 392         ( fact:fact*i,
 393           for j thru neqns do
 394           ( value:errcatch(diff(eqtn[j],x)),
 395             if value=[] then exitfor(order:i-1),
 396             eqtn[j]:first(value), value:errcatch(ratsimp(at(eqtn[j],xvar=xhat))),
 397             if value=[] then exitfor(order:i-1),
 398             ufun[j]:diff(ufun[j],x), value:first(value),
 399             atvalue(ufun[j], xvar=xhat, value),
 400             tfun[j]:tfun[j]+value*(x-xhat)^i/fact),
 401           if order#napprox then exitfor(false)),
 402         for i:neqns step -1 thru 1 do vlist:cons(tfun[i],vlist),
 403         printsoln(vlist, 'tailor, order, testsoln(vlist)))$
 404
 405 /* obtains expansion point for taylor series by solving b(x)=a(x).
 406   globals: xvar, xhat */
 407 getxhat(expr):=
 408         if mapatom(expr) then false else
 409         if opr(expr)#'?%integrate then for sub in expr do getxhat(sub) else
 410         if xhat='xhat then
 411           (xhat:mysolve(inpart(expr,3)-inpart(expr,4),xvar),
 412            if xhat=[] then throw(false) else xhat:first(xhat),
 413            if ratsubst(xhat, xvar, inpart(expr,3))#xhat then throw(false)) else
 414         if subst(xhat,xvar,expr)#0 then throw(false)$
 415
 416 /* fredholm-carleman series.   globals: napprox, kxu, ax, bx, fx, xvar, uvar */
 417 fredseries():=block([order, alpha, bta, top, bot, kind, tvar, eqtn, kxt],
 418         eqtn:first(eqnlist),
 419         alpha:catch(findint(eqtn)), if alpha=false then return(false),
 420         uvar:inpart(alpha,2), bta:lin(eqtn,alpha),
 421         if bta=false then return(false),
 422         kxu:lin(arg(alpha),subst(uvar,xvar,px)),
 423         if kxu=false then return(false),
 424         ax:inpart(alpha,3), bx:inpart(alpha,4), fx:bta[2],
 425         if kxu[2]#0 then fx:fx+myint(bta[1]*kxu[2],uvar,ax,bx),
 426         kxu:kxu[1]*bta[1], order:napprox, kind:'approximate,
 427         kxt:subst(tvar,uvar,kxu), bot:bta:1,
 428         top:alpha:rat(kxu,uvar,xvar),
 429         for i thru napprox do
 430         ( bta:-myint(ratsubst(uvar,xvar,alpha),uvar,ax,bx)/i,
 431           alpha:bta*kxu+myint(kxt*ratsubst(tvar,xvar,alpha),tvar,ax,bx),
 432           bot:bot+bta,
 433           if alpha=0 then exitfor(kind:[]),
 434           top:top+alpha),
 435         kxt:fx+myint(subst(uvar,xvar,fx)*top,uvar,ax,bx)/ratdisrep(bot),
 436         printsoln(kxt, 'fredseries, order, kind))$
 437
 438
 439 /* neumann series.   globals: eqnlist, pxlist, guesslist, napprox  */
 440 neumann():=block([order, newguess],
 441         apply('local, maplist('opr,pxlist)), order:napprox,
 442         if guesslist='none then guesslist:maplist('zeroint, eqnlist),
 443         for count thru napprox do
 444           (maplist('define, pxlist, guesslist),
 445            newguess:apply('ev,[eqnlist,'integrate]),
 446            if not freeof('?%integrate, newguess) then exitfor(order:count-1),
 447            guesslist:maplist('ratsimp,newguess)),
 448         printsoln(guesslist, 'neumann, order, testsoln(guesslist)))$
 449
 450
 451
 452 /* ---- the following functions apply only to 1st kind eqns */
 453
 454 /* fixed-limit finite-rank 1st kind.  globals: kxu, fx, uvar, xvar, ax, bx */
 455 flfrnk1st():=block([sf, cnt, intg, form, %c, unklist, eqlist,res],
 456
 457         local(%c), cnt:form:0, unklist:[],
 458         if not freeof(xvar,[ax,bx]) then return(false),
 459         intg:sumfactors(kxu,uvar),
 460         for term in intg do
 461           (cnt:cnt+1, form:form+term[2]*%c[cnt],
 462            unklist:cons(%c[cnt],unklist)),
 463         eqlist:[],
 464         sf:num(ratsimp(myint(kxu*form,uvar,ax,bx)-fx,xvar)),
 465         res:0, if opr(sf)#"+" then sf:[sf],
 466         for term in sf do
 467           if opr(term)="*" then
 468             (intg:partition(term,xvar),
 469              if intg[2]=1 then res:res+term else eqlist:cons(intg[1],eqlist))
 470           else if freeof(xvar,term) then res:res+term else error("error 1"),
 471         if res#0 then eqlist:cons(res,eqlist),
 472         if length(eqlist)#cnt then error("error 2"),
 473         solveandsubst(eqlist, unklist, subst(xvar,uvar,form), 'flfrnk1st))$
 474
 475 /* some=one*peace+two, where peace could be a product. */
 476 mydivide(some,peace,var):=block([res, one, two],
 477         one:two:0, if opr(some)#"+" then some:[some],
 478         for prt in some do
 479           (res:quotient(prt,peace),
 480            if res#0 and freeof(var,res) then one:one+res else two:two+prt),
 481         return([one,two]))$
 482
 483
 484 /* generalized abel method.  globals: kxu, ax, bx, fx, xvar, uvar */
 485 abel():=block([power,den,fun],
 486         if not freeof(xvar,ax) or bx#xvar or opr(kxu)#"^" then throw(false),
 487         power:-inpart(kxu,2), den:inpart(kxu,1),
 488         if not freeof(uvar, power) or opr(den)#"+" then throw(false),
 489         fun:partition(den,uvar),
 490         if not freeof(xvar, fun[2]) or
 491            ratsimp(subst(uvar,xvar,fun[1])+fun[2])#0 then throw(false),
 492         if sign(power)#'pos or sign(power-1)#'neg then throw(false),
 493         fun:myint(diff(-fun[2],uvar)*subst(uvar,xvar,fx)*den^(power-1),uvar,ax,bx),
 494         fun:sin(%pi*power)/%pi*diff(fun,xvar),
 495         printsoln(fun, 'abel, [], []))$
 496
 497
 498 /* firstkindseries.  globals: px, napprox, guess, eqn, px  */
 499 firstkindseries():=block([kind, order, correction, trial],
 500         apply('local,[opr(px)]), kind:'approximate, order:napprox,
 501         trial:if guess='none then zeroint(eqn) else guess,
 502         for i thru napprox do
 503           (apply('define,[px, trial]), correction:apply('ev,[eqn,'integrate]),
 504            if not freeof('?%integrate, correction) then exitfor(order:i-1),
 505            if (correction:ratsimp(correction))=0 then exitfor(kind:[]),
 506            trial:ratsimp(trial-correction)),
 507         printsoln(trial, 'firstkindseries, order, kind))$
 508
 509 /* kill(labels)$ */
 510 ttyoff:false$