1eval_when([translate,batch,demo,load,loadfile], 2dv(a)::=buildq([a],define_variable(a,'a,any)), 3dvl(a)::=buildq([a],define_variable(a,[],list)),ttyoff:true)$ 4 5/* (c) Copyright 1981 Massachusetts Institute of Technology */ 6 7/* dynamalloc:true$ */ 8define_variable(ieqnprint,true,boolean)$ 9define_variable(techlist,'[first,all,vlfrnk,transform,collocate, 10 flfrnk2nd,tailor,fredseries,neumann, 11 flfrnk1st,abel,firstkindseries],list)$ 12dv(tech)$ 13dv(px)$ 14dv(eqn)$ 15dv(fx)$ 16dv(uvar)$ 17dv(xvar)$ 18dv(kxu)$ 19dv(ax)$ 20dv(bx)$ 21dv(iesoln)$ 22dv(napprox)$ 23dv(pxlist)$ 24dv(eqnlist)$ 25dv(inteqn)$ 26dv(lowlim)$ 27dv(pofx)$ 28dv(iclist)$ 29dvl(unklist)$ 30dv(highlim)$ 31dv(%c)$ 32dv(eqno)$ 33dvl(ueqnlist)$ 34dv(x)$ 35dv(xhat)$ 36dv(opr)$ 37dv(guesslist)$ 38dv(guess)$ 39 40ieqn([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 47opr(e):=if mapatom(e) then 'none else inpart(e,0)$ 48arg(e):=inpart(e,1)$ 49alias(exitfor, return)$ 50 51/* identifies functions unknown to macsyma */ 52unfun(e):=is(?getchar(e,1)='?\$)$ 53 54/* ***** main program *****/ 55ieqn1(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 */ 90try(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 */ 99firstkind():=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 */ 116findint(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 */ 122printsoln(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 134alias(myint, integrate)$ 135 136myprint(msg1,msg2):=if ieqnprint then print(msg1,msg2)$ 137 138/* if expr = a*var+b returns [a,b] else false */ 139lin(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 */ 144mysolve(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. */ 151solvesys(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 */ 157sumfactors(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 */ 171frnk(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 */ 193plussplit(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 */ 200solveandsubst(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 */ 212testsoln(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 */ 222mypartition(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 */ 237zeroint(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 */ 246vlfrnk():=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 */ 276getics():=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 284be xvar and lower limit must be a constant. 285 globals: initial, lowlim, unklist, xvar */ 286vconvert(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 */ 304transform():=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 */ 316collocate():=block([lowlim, highlim, unklist, %c, elist, form, incr, point, 317 name, listeqns, solnlist], 318 apply('local, cons('%c,maplist('opr,pxlist))), lowlim:'minf, 319highlim:'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 */ 337getlimits(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) */ 346approx(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 */ 351flfrnk2nd():=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 occurences 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 (indepedent variable) */ 368fconvert(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 */ 383tailor():=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 */ 407getxhat(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 */ 417fredseries():=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 */ 440neumann():=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 */ 455flfrnk1st():=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. */ 476mydivide(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 */ 485abel():=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 */ 499firstkindseries():=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)$ */ 510ttyoff:false$ 511