1%********************************************************************* 2module integration$ 3%********************************************************************* 4% Routines for integration of pde's 5% Authors: Andreas Brand 1993 1995 6% Thomas Wolf since 1993 7 8% BSDlicense: ***************************************************************** 9% * 10% Redistribution and use in source and binary forms, with or without * 11% modification, are permitted provided that the following conditions are met: * 12% * 13% * Redistributions of source code must retain the relevant copyright * 14% notice, this list of conditions and the following disclaimer. * 15% * Redistributions in binary form must reproduce the above copyright * 16% notice, this list of conditions and the following disclaimer in the * 17% documentation and/or other materials provided with the distribution. * 18% * 19% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * 20% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * 21% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * 22% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR CONTRIBUTORS BE * 23% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR * 24% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF * 25% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS * 26% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN * 27% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * 28% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * 29% POSSIBILITY OF SUCH DAMAGE. * 30%****************************************************************************** 31 32symbolic procedure ldlist(p,f,vl)$ 33% provides a reverse list of leading derivatives of f in p, vl is list 34% of variables 35begin scalar a$ 36 if not atom p then 37 if member(car p,list('expt,'plus,'minus,'times,'quotient,'df,'equal)) 38 then << 39 if (car p='plus) or (car p='times) or 40 (car p='quotient) or (car p='equal) then 41 <<p:=cdr p$ 42 while p do 43 <<a:=diffincl(a,ldlist(car p,f,vl))$ 44 p:=cdr p>> 45 >> else 46 if car p='minus then a:=ldlist(cadr p,f,vl) else 47 if car p='expt then % if numberp caddr p then 48 a:=ldlist(cadr p,f,vl) else % fuehrende Abl. aus der Basis 49 % else a:=nil 50 if car p='df then if cadr p=f then 51 <<p:=cddr p; 52 while vl do 53 <<a:=cons(dfdeg(p,car vl),a); 54 vl:=cdr vl>>; 55 a:=list a 56 >> 57 >>$ 58 return a 59end$ 60 61symbolic procedure diffincl(a,b)$ 62% a,b are lists of leading derivatives which are to be united 63begin 64 scalar p; 65 while a and b do 66 <<a:=ddroplow(a,car b); 67 if car a then p:=cons(car a,p); 68 a:=cdr a; 69 b:=cdr b>>; 70 return 71 if null a then if p then nconc(p,b) 72 else b 73 else if p then a:=nconc(p,a) 74 else a 75end$ 76 77symbolic procedure ddroplow(a,cb)$ 78% loescht Elemente von a, von denen cb eine Ableitung ist, loescht cb, 79% wenn ein Element von a eine Ableitung von cb ist 80begin 81 scalar h; 82 return 83 if null a then list(cb) 84 else 85 <<h:=compdiffer(car a,cb); 86 if numberp(h) then if h>0 then cons(nil,a) % car a=df(cb,.. 87 else ddroplow(cdr a,cb) % cb=df(car a,.. 88 else <<h:=ddroplow(cdr a,cb); % neither 89 cons(car h,cons(car a,cdr h))>> 90 >> 91end$ 92 93symbolic procedure compdiffer(p,q)$ 94% finds whether p is a derivative of q or q of p or neither 95begin 96 scalar p!>q,q!>p; 97 while p and ((null p!>q) or (null q!>p)) do 98 <<if car p>car q then p!>q:=t else % compare orders of derivatives 99 if car p<car q then q!>p:=t; 100 p:=cdr p;q:=cdr q 101 >>; 102 return 103 if q!>p then if p!>q then nil % neither 104 else 0 % q is derivative of p 105 else if p!>q then 2 % p is derivative of q 106 else 1 % p equal q 107end$ 108 109 110symbolic procedure ldintersec(a)$ 111% determines the intersection of derivatives in the list a 112begin 113 scalar b; 114 return 115 if null a then nil else 116 <<b:=car a;a:=cdr a; 117 while a do 118 <<b:=isec(b,car a)$ 119 a:=cdr a 120 >>; 121 b 122 >> 123end$ 124 125 126symbolic procedure isec(ca,b)$ 127% determines the minimum derivatives between ca and b 128begin 129 scalar newb; 130 while ca do 131 <<newb:=cons(if car b<car ca then car b else car ca, newb); 132 ca:=cdr ca;b:=cdr b 133 >>; 134 return reverse newb 135end$ 136 137 138symbolic procedure disjun(a,b)$ 139<<while a do 140 if (car a neq 0) and (car b neq 0) then a:=nil 141 else <<a:=cdr a;b:=cdr b>>; 142 if b then nil else t 143>>$ 144 145 146symbolic procedure ddrophigh(a,cb)$ 147% loescht Elemente von a, die Ableitung von cb sind, 148% loescht cb, wenn cb Ableitung eines Elements von a ist oder =a ist, 149% haengt cb ansonsten an 150begin 151 scalar h; 152 return 153 if null a then list(cb) 154 else 155 <<h:=compdiffer(car a,cb); 156 if numberp(h) then if h<2 then a % cb is deriv or = car a 157 else ddrophigh(cdr a,cb) % car a=df(cb,.. 158 else cons(car a,ddrophigh(cdr a,cb)) % neither 159 >> 160end$ 161 162 163symbolic procedure elibar(l1,l2,lds)$ 164begin 165 scalar found1,found2,h; 166 % found1=t if an LD=l1 is found, found2=t if contradiction found 167 while lds and (not found2) do 168 <<if car lds=l1 then found1:=t else 169 <<h:=compdiffer(l2,car lds); 170 if (null h) or (h = 2) then found2:=t 171 >>; 172 lds:=cdr lds 173 >>; 174 return found1 and (not found2) 175end$ 176 177symbolic procedure intminderiv(p,ftem,vlrev,maxvanz,nfsub)$ 178% yields a list {nv1,nv2, ... } such that nvi is the minimum 179% of the highest derivatives w.r.t. vi of all the leading derivatives 180% of all functions of ftem which are functions of all maxvanz variables. 181% Only those are kept for which nvi>0. 182% further a list ((f1 ld's of f1) (f2 ld's of f2)...), 183begin scalar l,a,listoflds$ 184 while ftem do 185 <<if (maxvanz=(fctlength car ftem)) or (nfsub=0) then 186 <<l:=ldlist(p,car ftem,vlrev); 187 listoflds:=cons(cons(car ftem,l),listoflds)$ 188 a:=if a then ldintersec(cons(a,l)) 189 else ldintersec(l) 190 >>$ 191 ftem:=cdr ftem 192 >>$ 193 return list(a,listoflds) 194end$ 195 196 197symbolic procedure potintegrable(listoflds)$ 198begin 199 scalar h,l1,l2,f,n,lds,f1,f2$ 200 if tr_genint then write "Does a potential exist?"$ 201 %------- Determining 2 minimal sets of integration variables 202 % There must be two disjunct LDs such that all others are in their 203 % ideal. This is necessary if we want to eliminate 2 functions. 204 h:=listoflds; 205 l1:=nil; 206 while h do 207 <<l2:=cdar h; % the list of LDs for the function caar h 208 while l2 do 209 <<l1:=ddrophigh(l1,car l2)$ 210 l2:=cdr l2>>; 211 h:=cdr h 212 >>; 213 return 214 if length l1 neq 2 then nil else 215 if not disjun(car l1,cadr l1) then nil else 216 % if there would be an overlap between l1 and l2 then it would have 217 % to be integrated for elimination but it cannot --> no overlap 218 % possible 219 <<% selecting interesting functions for which one LD is = l1 and all 220 % others are derivatives of l2 or equal l2 and vice versa. Two 221 % necessary (one with an LD=l1 and one with an LD=l2) from which at 222 % least one function (f) has no further LD. 223 % Exception: 2 functions with each 2 LDs equal to (l1,l2) (these 224 % functions are counted by n) 225 l2:=cadr l1;l1:=car l1; 226 f:=nil;f1:=nil;f2:=nil;n:=0; 227 h:=listoflds; 228 while h and ((not f1) or (not f2) or ((not f) and (n neq 2))) do 229 <<lds:=cdar h; 230 if (not f1) or (not f) then 231 if elibar(l1,l2,lds) then 232 <<f1:=cons(caar h,f1); 233 if length lds = 1 then f:=caar h else 234 if length lds = 2 then 235 if (car lds = l2) or (cadr lds = l2) then n:=n+1 236 >>; 237 if (not f2) or (not f) then 238 if elibar(l2,l1,lds) then 239 <<f2:=cons(caar h,f2); 240 if length lds = 1 then f:=caar h 241 >>; 242 h:=cdr h 243 >>; 244 if f1 and ((n>1) or (f2 and f)) then list(l1,l2) 245 else nil 246 >> 247end$ % of potintegrable 248 249symbolic procedure vlofintlist(vl,intlist)$ 250% provides a list of elements of vl for which the corresponding 251% elements of intlist are non-zero 252begin scalar a; 253 while intlist do 254 <<if (car intlist) and (not zerop car intlist) then a:=cons(car vl,a); 255 vl:=cdr vl; 256 intlist:=cdr intlist 257 >>; 258 return a 259end$ 260 261symbolic procedure vlofintfaclist(vl,intfaclist)$ 262% determining the list of all variables of vl in intfaclist 263begin scalar e1,a; 264 for each e1 in vl do 265 if not my_freeof(intfaclist,e1) then a:=cons(e1,a); 266 return a 267end$ 268 269symbolic procedure multipleint(intvar,ftem,q,vari,vl,genflag, 270 potflag,partial,doneintvar)$ 271% multiple integration of q wrt. variables in vl, max. number of 272% integrations specified by intvar 273% integrating factors must not depend on doneintvar, doneintvar is 274% a list of integration variables so far 275% partial=t then as much as possible of an expression is integrated 276% even if there is a remainder 277begin 278 scalar pri,vlcop,v,nmax,geni,intlist,iflag,n,nges,newcond, 279 intfaclist,ph,pih,qh,qih,intfacdepnew,intfacdep$ 280 % intfacdep is a list of variables on which factors of integration 281 % depend so far, other than the integration variable in their 282 % integration --> no integration wrt. these variables by potint 283 % because there the diff. operators wrt. to different variables 284 % need not commute because the integrations are not done 285 286 % pri:=t$ 287 if (not vari) and (zerop q) then return nil; 288 nges:=0; 289 vlcop:=vl; 290 pih:=t; 291 292 % Splitting of the equation into the homogeneous and inhomogeneous 293 % part which is of advantage for finding integrating factors 294 q:=splitinhom(q,ftem,vl)$ 295 qh:=car q; qih:=cdr q; q:=nil; 296 297 while (vari or vlcop) and (pih or (not potflag)) do 298 %------- if for potflag=t one variable can not be integrated the 299 %------- maximal number of times (nmax) then immediate stop because 300 %------- then no elimination of two functions will be possible 301 << %-- The next integration variable: v, no of integrations: nmax 302 if vari then <<v:=vari;nmax:=1>> 303 else <<v:=car vlcop; vlcop:=cdr vlcop; 304 nmax:=car intvar; intvar:=cdr intvar>>; 305 306 if zerop nmax then intlist:=cons(nil,intlist) 307 else 308 <<if pri then write"anf: intvar=",intvar," vari=",vari," q=",q, 309 " v =",v," vl =",vl$ 310 if vari and (not member(v,vl)) then 311 <<%qh :=reval list('INT,qh,v)$ 312 qh :=err_catch_int(qh,v)$ % changed 23.1.08 313 if null qh then iflag:=nil else % changed 23.1.08 314 if freeof(qh,'int) then << 315 %qih:=reval list('INT,qih,v)$ 316 qih:=err_catch_int(qih,v)$ % changed 23.1.08 317 iflag:=if null qih then nil else % changed 23.1.08 318 if freeint_ and 319 (null freeof(qih,'int)) then nil else 320 if freeabs_ and 321 (null freeof(qih,'abs)) then nil else << 322 intlist:=cons(list(1),intlist)$ 323 'success>>$ 324 if pri then <<write"232323 qh=",qh;terpri(); 325 write"qih=",qih;terpri()>> 326 >> 327 >> else 328 <<n:=0$ 329 if pri then write"333"$ 330 intfaclist:=nil; %-- the list of integr. factors in v-integr. 331 if potflag or my_freeof(intfacdep,v) then 332 % otherwise v-integration not allowed because one integrating 333 % factor already depends on v 334 % for potflag=t this `commutativity'-demand plays no role 335 repeat << %--- max nmax integrations of qh and qih wrt. v 336 if pri then <<write"444 vor intpde:"$eqprint q$terpri()$ 337 write"potflag=",potflag," v=",v, 338 " ftem=",ftem>>$ 339 % At first trying a direct integration of the homog. part qh 340 ph:=intpde(qh,ftem,vl,v,partial)$ % faster if potflag=nil 341 if pri then <<write"nach intpde(qh):"$deprint ph>>$ 342 343 %------ At first the integration of the homogeneous part 344 intfacdepnew:=intfacdep; 345 if ph and (partial or (zerop cadr ph)) then << 346 %---- For the homogen. part cadr ph must be zero 347 intfaclist:=cons(1,intfaclist); 348 ph:=car ph; 349 if pri then <<write"565656 ph=",ph;terpri()>>; 350 >> else 351 if vari then ph:=nil 352 else 353 if facint_ then << 354 ph:=findintfac(list(qh),ftem,vl,v,doneintvar,intfacdep, 355 not zerop n,not potflag); 356 % factorize before ivestig., no report of int. factors 357 if ph then << %--- Complete integr. of qh was possible 358 if pri then <<write"of the homogeneous part"$terpri()>>$ 359 %--- update the list of variables on which all integr. 360 %--- factors depend apart from the integration variable 361 intfacdepnew:=caddr ph; 362 %--- extend the list of integrating factors, cadr ph 363 %--- is a list of integr. factors, here only one 364 intfaclist:=cons(caadr ph,intfaclist); 365 %--- multiply the inhomogeneous part with integ. factor 366 qih:=reval reval reval list('times,car intfaclist,qih); 367 if pri then <<write"454545 qih=",qih;terpri()>>; 368 ph:=car ph % the integral of qh 369 >> 370 >>; 371 372 %------ Now the integration of the inhomogeneous part 373 if not ph then pih:=nil %--- no integration possible 374 else << 375 if zerop qih then pih:=list(0,0) else 376 pih:=intpde(qih,ftem,vl,v,t)$ % partial set =t on 14.6.04 377 % to generalize integrate for inhom. terms 378 if print_ and null pih then << 379 terpri()$ 380 write"Inhomogeneous part: "$ 381% mathprint qih$ %@!@!@!@! 382 type_pre_ex qih$ 383 write"can not be integrated explicitly wrt. ",v$ 384 >>$ 385 if pri then <<write"nach intpde(qih):",pih$terpri()$ 386 write"genflag=",genflag$terpri()>>$ 387 if pih then 388 if zerop cadr pih then 389 <<qih:=car pih$n:=add1 n$iflag:='success$ 390 if pri then write" success "$ 391 >> 392 else if not genflag then pih:=nil 393 else 394 <<if pri then write"555"$ 395 geni:=partint(cadr pih,smemberl(ftem,cadr pih), 396 vl,v,genflag)$ 397 if geni then 398 <<qih:=reval list('plus,car pih,car geni)$ 399 n:=add1 n$ 400 ftem:=union(fnew_,ftem)$ 401 newcond:=append(cdr geni,newcond)$ % additional de's 402 if pri then 403 <<terpri()$write"nach gen newcond:",newcond>>$ 404 iflag:='genint 405 >> else 406 if partial then qih:=car pih 407 else pih:=nil 408 >>; 409 if pih then << 410 if pri then write"AAA"$ 411 qh:=ph; 412 if (not potflag) and (n neq 1) then 413 intfacdep:=intfacdepnew 414 %-The first integr. factor of all v-integrations does not 415 % depend on any earlier integration variables and can 416 % therefore be taken out of all integrals --> no incease 417 % of intfacdep for n=1. 418 %-For potential integration the integration variables and 419 % extra-dependencies-variables of integr. factors need not 420 % be disjunct therefore the intfacdep-update only for 421 % not-potflag 422 >> else << 423 if pri then write"BBB"$ 424 % inhomogeneous part can not be integrated therefore 425 % reversing the succesful integration of the hom. part 426 if car intfaclist neq 1 then 427 qih:=reval list('quotient,qih,car intfaclist); 428 intfaclist:=cdr intfaclist 429 >>; 430 >>; %-- end of the integration of the inhomog. part 431 if pri then write"n=",n," nmax=",nmax," not pih=",not pih$ 432 >> until (n=nmax) or (not pih)$ %---- end of v-integration 433 %------- variables of done integrations are collected for 434 %------- determining integrating factors in later integr. 435 if not zerop n then doneintvar:=cons(v,doneintvar)$ 436 nges:=nges+n; 437 intlist:=cons(intfaclist,intlist) 438 >>$ % of not ( vari and (not member(v,vl))) 439 if vari then <<vari:=nil;vlcop:=nil>>; 440 if pri then write"ende: intvar=",intvar," vari=",vari, 441 " qh=",qh," qih=",qih$ 442 >> % of (nmax neq zero) 443 >>$ % of ( while (vari or vlcop) and (pih or (not potflag)) ) 444 445 % intlist and its elements intfaclist are in the inverse order 446 % to vl and the sequence of integrations done 447 q:=reval list('plus,qh,qih)$ %--- adding homog. and inhomog. part 448 if pri then <<terpri()$write"\\\ newcond:"$deprint newcond; 449 write"multiple result:",if null iflag then nil 450 else list(q,intlist,newcond,nges) 451 >>; 452 return if null iflag then nil 453 else list(q,intlist,newcond,nges) 454end$ % of multipleint 455 456symbolic procedure uplistoflds(intlist,listoflds)$ 457begin 458 scalar f,h1,h2,h3,h4,lds,itl; 459 while listoflds do 460 <<f:=caar listoflds; 461 lds:=cdar listoflds; 462 listoflds:=cdr listoflds; 463 h2:=nil; % h2 becomes the new list of lds of f 464 while lds do 465 <<h3:=car lds; lds:=cdr lds; 466 itl:=intlist; 467 h4:=nil; % h4 becomes one new ld of f 468 while h3 do 469 <<h4:=cons(car h3 - if null car itl then 0 470 else length car itl, h4); 471 h3:=cdr h3;itl:=cdr itl 472 >>; 473 h2:=cons(reverse h4,h2) 474 >>; 475 h1:=cons(cons(f,h2),h1) 476 >>; 477 return h1 % updated listoflds 478end$ % of uplistoflds 479 480symbolic procedure ProportionalityConditions(ex,fl,x)$ 481% This procedure collects cases that lead to at least 482% two coefficients of two elements of fl in ex being 483% proportional to each other by an x-independent multiplier. 484if fl then 485begin scalar flo,f1,f2,flc,s1,s2,c1,condi$ 486 flo:=fl$ 487 while cdr fl do << 488 f1:=car fl; fl:=cdr fl; 489 s1:=coeffn(ex,f1,1)$ 490 if not zerop s1 then << 491 flc:=fl; 492 while flc do << 493 f2:=car flc; flc:=cdr flc; 494 s2:=coeffn(ex,f2,1)$ 495 if not zerop s2 then << 496 497 % condition for checking a special case: 498 c1:=reval {'df,{'quotient,s1,s2},x}; 499 c1:=simplifySQ(simp c1,ftem_,t,nil,t)$ 500 if (c1 neq {(1 . 1)}) and (not freeoflist(c1,ftem_)) then 501 condi:=union(c1,condi) 502 >> 503 >> 504 >> 505 >>$ 506 return condi 507end$ 508 509symbolic procedure addintco(q, ftem, ifac, vl, vari)$ 510begin scalar v,f,l,vl1,j,ftemcp,fnewcp; 511 % multi.ing factors to the constants/functions of integration 512 ftemcp:=ftem; 513 if zerop q then l:=1 514 else 515 <<ftem:=fctsort ftem$ 516 while ftem do 517 if fctlength car ftem<length vl then ftem:=nil 518 else if fctlinear(q,f) then 519 <<f:=car ftem$ ftem:=nil>> else 520 ftem:=cdr ftem$ 521 if f then 522 <<l:=lderiv(q,f,fctargs f)$ 523 l:=reval coeffn(q,reval car l,cdr l)$ 524 % l is a coeffient of the leading derivative. By multiplying to the 525 % constant(function) of integration, the factor will dissappear when 526 % a substitution will be made. This may be dangerous because it 527 % may hide the case when the factor becomes zero (although no division 528 % may be performed through a factor which might become zero). But even 529 % when no division is performed then the constant of integration may 530 % dissappear when this factor becomes zero and thus solutions be lost. 531 % Therefore: 532 if not freeoflist(l,ftem) then l:=1 533 >> else l:=1 534 >>; 535 % the constants and functions of integration 536 if vari then q:=list('plus,q,intconst(l,vl,vari,list(1))) 537 % The coefficient is 1 so no testing of case splitting due to 538 % specific parameter values is done. 539 else 540 <<vl1:=vl; 541 while vl1 and null j do % j = list of case distinctions to do 542 <<v:=car vl1; vl1:=cdr vl1; 543 fnewcp:=fnew_; 544 if car ifac then 545 q:=list('plus,q,intconst(l,vl,v,car ifac))$ 546 % l..product of factors in the coefficient of the function to be 547 % eliminated, car ifac .. list of integrating factors 548 549 % All integrations wrt. v are done and now the independence of special 550 % solutions and the absence of singularities has to be tested. Both 551 % would give rise to more case distinctions. 552 % singularity check: 553 j:=zero_den(q,ftemcp); 554 % proportionality check: 555 if null j then j:=ProportionalityConditions(q,setdiff(fnew_,fnewcp),v)$ 556 557 ifac:=cdr ifac; 558 >> 559 >>$ 560 561 return 562 if null j then reval q 563 else << % case distinctions need to be made 564 if freeoflist(j,fnew_) then 565 % Otherwise conditions contain new constants of integration which 566 % are not known later after the integration module is left 567 % (see 30.1.2015) 568 for each h in j do 569 to_do_list:=cons(list('split_into_cases,h),to_do_list); 570 nil 571 >> 572 573end$ % of addintco 574 575symbolic procedure integratepde(p,ftem,vari,genflag,potflag)$ 576% Generalized integration of the expression p 577% if not genflag then "normal integration" 578% Equation p must not be directly separable, otherwise the depen- 579% dencies of functions of integration on their variables is wrong, 580% i.e. no dependence on explicit variables 581% ftem are all functions from the `global' ftem which occur in p, i.e. 582% ftem:=smemberl(ftem,p)$ 583% if vari=nil then integration w.r.t. all possible variables within 584% the equation 585% else only w.r.t. vari one time 586 587begin 588 scalar vl,vlrev,v,intlist, 589 ili1a,ili2a,maxvanz,fsub,h,hh,nfsub,iflag,newcond, 590 n1,n2,pot1,pot2,p1,p2,listoflds,secnd,ifac0, 591 ifac1a,ifac1b,ifac2a,ifac2b,cop,v1a,v2a,pri,aic,pnew$ 592 593% pri:=t; 594 if pri then <<terpri()$write"Start Integratepde">>$ 595 vl:=argset ftem$ 596 vlrev:=reverse vl; 597 if vari then <<potflag:=nil; 598 if zerop p then iflag:='success>> 599 else 600 <<%------- determining fsub=list of functions of all variables 601 maxvanz:=length vl$ 602 fsub:=nil; 603 h:=ftem; 604 while h do 605 <<if fctlength car h=maxvanz then 606 fsub:=cons(car h,fsub)$ 607 h:=cdr h 608 >>$ 609 nfsub:=length fsub$ % must be >1 for potential-integration 610 h:=intminderiv(p,ftem,vlrev,maxvanz,nfsub)$ % fsub is also for below 611 intlist:=car h$ 612 %-- list of necessary integrations of the whole equation to solve 613 %-- for a function of all variables 614 listoflds:=cadr h$ %-- list of leading derivatives 615 >>$ 616 if pri then <<terpri()$ 617 write"complete integrations:",intlist," for:",vl>>; 618 %-- n1 is the number of integrations which must be done to try 619 %-- potential integration which must enable to eliminate 2 functions 620 %-- n2 is the number of integrations actually done 621 n1:=for each h in intlist sum h; 622 if (not vari) and (zerop n1) then 623 <<n2:=0; 624 if potflag then % else not necessary 625 for h:=1:(length vl) do ifac0:=cons(nil,ifac0) 626 >> else 627 <<if tr_genint then 628 <<terpri()$write "integration of the expression : "$ 629 eqprint p>>$ 630 if pri then 631 <<terpri()$write"at first all multiple complete integration">>; 632 %-- At first if possible n2 integrations of the whole equation 633 h:=multipleint(intlist,ftem,p,vari,vl,genflag,nil,nil,nil)$ 634 % potflag=nil, partial=nil, doneintvar=nil 635 if h then 636 <<p:=car h; 637 ifac0:=cadr h; % ifac0 is the list of lists of integr. factors 638 newcond:=caddr h; 639 n2:=cadddr h; % number of done integrations 640 % doneintvar & intfacdep for the two halfs of integrations 641 % from the two parts of ifac0 642 h:=nil; 643 iflag:='success; 644 >> else n2:=0; 645 ftem:=union(fnew_,ftem)$ 646 >>; 647 %------------ Existence of a potential ? 648 if (n1=n2) and potflag and (nfsub>1) then 649 %---- at least 2 functions to solve for 650 <<if not zerop n2 then %---- update listoflds 651 listoflds:=uplistoflds(reverse ifac0,listoflds)$ 652 if pri then <<terpri()$write"uplistoflds:",listoflds>>$ 653 if h:=potintegrable(listoflds) then 654 <<ili1a:=car h; ili2a:=cadr h; 655 % The necess. differentiations of the potential 656 if pri then 657 <<terpri()$write"potintegrable:",ili1a," ",ili2a>>$ 658 659 if pri then <<write"+++ intlist=",intlist, 660 " ili1a=",ili1a, 661 " ili2a=",ili2a>>$ 662 %-- distributing the integrating factors of ifac0 among 663 %-- the two lists ifac1b and ifac2b which are so far nil 664 %-- such that (ifac1b and ili1a are disjunct) and 665 %-- (ifac2b and ili2a are disjunct) 666 v1a:=vlofintlist(vl,ili1a); 667 v2a:=vlofintlist(vl,ili2a); 668 669 hh:=t; 670 cop:=reverse ifac0; 671 ifac1a:=ili1a;ifac2a:=ili2a; 672 while hh and cop do << 673 % cop is a list of lists of integr. factors 674 if car cop then h:=vlofintfaclist(vl,cdar cop) 675 else h:=nil; 676 if freeoflist(h,v2a) and (car ifac2a=0) then << 677 ifac1b:=cons( nil, ifac1b); 678 ifac2b:=cons( reverse car cop, ifac2b) 679 >> else 680 if freeoflist(h,v1a) and (car ifac1a=0) then << 681 ifac2b:=cons( nil, ifac2b); 682 ifac1b:=cons( reverse car cop, ifac1b) 683 >> else 684 if car cop then hh:=nil; 685 ifac1a:=cdr ifac1a; 686 ifac2a:=cdr ifac2a; 687 cop:=cdr cop; 688 >>; 689 % the elements of ifac1b,ifac2b are in reverse order to 690 % ifac1a,ifac2a and are in the same order as vl, also 691 % the elements in the infac-lists are in inverse order, 692 % i.e. in the order the integrations have been done 693 if pri then <<terpri()$ 694 write "ifac1a=",ifac1a," ifac1b=",ifac1b, 695 " ifac2a=",ifac2a," ifac2b=",ifac2b >>$ 696 697 %-- lists of integrations to be done to both parts 698 if hh then 699 repeat % possibly a second try with part2 integrated first 700 <<n1:=for each n1 in ili1a sum n1; 701 % n1 .. number of remaining integrations of the first half 702 p1:=multipleint(ili1a,ftem,p,nil,vl,genflag,t,t, 703 % potflag=t, partial=t 704 union(vlofintlist(vl,ili2a), 705 vlofintlist(vl,ifac1b)))$ 706 % that the variables of integration are not in ifac1b 707 % was already checked. Only restriction: the integrating 708 % factors must not depend on the last argument. 709 710 ftem:=union(fnew_,ftem)$ 711 if p1 then << 712 ifac1a:=cadr p1; 713 % ifac1a is now the list of integrating factors 714 if newcond then newcond:=nconc(newcond,caddr p1) 715 else newcond:=caddr p1; 716 if pri then <<terpri()$write"mul2: newcond=",newcond>>$ 717 n2:=cadddr p1; 718 p1:=car p1 719 >>; 720 if p1 and (n1=n2) then 721 %--- if the first half has been integrated suff. often 722 <<%--- integrating the second half sufficiently often 723 n1:=for each n1 in ili2a sum n1; 724 % calculation of the 2. part which is not contained in p1 725 p2:=p1; 726 cop:=ifac1a; hh:=vlrev; % because ifac1a is reversed 727 while cop do << 728 h:=car cop;cop:=cdr cop; 729 v:=car hh;hh:=cdr hh; 730 % h is the list of integrating factors of the v-integr. 731 while h do << 732 p2:=reval list('quotient,list('df,p2,v),car h); 733 h:=cdr h 734 >> 735 >>; 736 p2:=reval reval list('plus,p,list('minus,p2)); 737 p2:=multipleint(ili2a,ftem,p2,nil,vl,genflag,t,nil, 738 % potflag=t, partial=nil 739 union(vlofintlist(vl,ili1a), 740 vlofintlist(vl,ifac2b)))$ 741 ftem:=union(fnew_,ftem)$ 742 if p2 then << 743 ifac2a:=cadr p2; 744 % ifac2a is now list of integrating factors 745 if newcond then newcond:=nconc(newcond,caddr p2) 746 else newcond:=caddr p2; 747 if pri then <<terpri()$write"mul3: newcond=",newcond>>$ 748 n2:=cadddr p2; 749 p2:=car p2 750 >>; 751 if p2 and (n1=n2) then 752 % if the second half has been integrated sufficiently often 753 <<% can both halfes be solved for different functions 754 % i.e. are they algebraic and linear in different functions? 755 pot1:=nil; 756 pot2:=nil; 757 for each h in fsub do << 758 if ld_deriv_search(p1,h,vl) = (nil . 1) then 759 pot1:=cons(h,pot1); 760 if ld_deriv_search(p2,h,vl) = (nil . 1) then 761 pot2:=cons(h,pot2); 762 >>$ 763 if (null not_included(pot1,pot2)) or 764 (null not_included(pot2,pot1)) then p2:=nil 765 >>$ 766 if p2 and (n1=n2) then 767 <<% ifac1b,ifac2b are in reverse order to ifac1a,ifac2a! 768 pot1:=newfct(fname_,vl,nfct_)$ % the new potential fct. 769 pot2:=pot1; 770 nfct_:=add1 nfct_$ 771 fnew_:=cons(pot1,fnew_); 772 flin_:=fctinsert(pot1,flin_)$ 773 v:=vlrev; 774 while v do 775 <<cop:=car ifac1a; ifac1a:=cdr ifac1a; 776 while cop do << 777 pot1:=reval list('quotient,list('df,pot1,car v),car cop); 778 cop:=cdr cop 779 >>; 780 cop:=car ifac2a; ifac2a:=cdr ifac2a; 781 while cop do << 782 pot2:=reval list('quotient,list('df,pot2,car v),car cop); 783 cop:=cdr cop 784 >>; 785 v:=cdr v; 786 >>; 787 pnew:=addintco(list('plus,p1,reval pot2), 788 ftem,ifac1b,vlrev,nil)$ 789 % This value is called pnew and not p because if secnd=nil then 790 % if pnew=nil or aic=nil below then another try can be done 791 % by integrating in a different order and then the former value 792 % of p is still needed. 793 % BUT, if pnew=nil or aic=nil then case-distinctions have 794 % already been booked in addintco() (in to_do_list) and then 795 % repeating in a different order would perhaps book the same 796 % case distinctions again which would not be good --> we set 797 % secnd:=t. 798 799 if null pnew then secnd:=t 800 else << 801 aic:=addintco(list('plus,p2,list('minus,reval pot1)), 802 ftem,ifac2b,vlrev,nil)$ 803 if null aic then secnd:=t 804 else << 805 newcond:=cons(aic,newcond)$ % vari=nil 806 iflag:='potint % i.e. success by integration with potential 807 >> 808 >> 809 >> 810 ;if pri then write":::"$ 811 >>; 812 % Before the following assignment it is secnd=t if this is the 813 % second time this loop is run and then no more run. 814 secnd:=not secnd; 815 % retry in different order of integration, p is still the same 816 if (iflag neq 'potint) and secnd then 817 <<cop:=ili1a;ili1a:=ili2a;ili2a:=cop>> 818 >> until (iflag eq 'potint) % success 819 or (not secnd) % no success (iflag=nil) 820 >>; 821 >>$ 822 823 %--------- returning the result 824 return if null iflag then nil 825 else 826 <<if iflag='potint then % pnew with contants of integration is computed. 827 else % add constants of integration 828 pnew:=addintco(p, ftem, % the completely reversed ifac0 829 <<h:=nil; 830 while ifac0 do <<h:=cons(reverse car ifac0,h); 831 ifac0:=cdr ifac0>>; 832 h 833 >>, vl, vari)$ 834 835 % If the terms involving constants of integration are not unique 836 % because their structure depends on the value of parameters 837 % (like int(x^n,x) = x^(n+1)/(n+1) or log(x) for n=-1) then 838 % case distinctions have been booked and then p=nil: 839 if null pnew then nil 840 else << 841 if pri then <<terpri()$write"ENDE INTEGRATEPDE"$ 842 deprint(cons(pnew,newcond))>>$ 843 cons(pnew,newcond) 844 >> 845 >> 846end$ % of integratepde 847 848 849symbolic procedure intpde(p,ftem,vl,x,potint)$ 850begin scalar ft,ip,h,itgl1,itgl2$ 851 852 if potint or null lin_problem then return intpde_(p,ftem,vl,x,potint)$ 853 % test null lin_problem added 28.5.03 854 855 % ft are functions of x 856 for each h in ftem do 857 if not freeof(assoc(h,depl!*),x) then ft:=cons(h,ft); 858 859 ip:=int_partition(p,ft,x)$ 860 if null cdr ip then return intpde_(p,ftem,vl,x,potint)$ 861 862 while ip do << 863 h:=intpde_(car ip,ftem,vl,x,potint)$ 864 if null h then << 865 ip:=nil; 866 itgl1:=nil; 867 itgl2:=nil 868 >> else << 869 % itgl1:=cons(car h,itgl1); 870 % itgl2:=cons(cadr h,itgl2); 871 itgl1:=nconc(list(car h),itgl1); 872 itgl2:=nconc(list(cadr h),itgl2); 873 ip:=cdr ip 874 >>; 875 >>$ 876 877 if itgl1 then <<itgl1:=reval cons('plus,itgl1); 878 itgl2:=reval cons('plus,itgl2) >>$ 879 880 return if null itgl1 then nil 881 else {itgl1,itgl2} 882end$ 883 884 885symbolic procedure drop_x_dif(der,x)$ 886begin scalar dv,newdv$ 887 % der is a derivative like {'df,f,u,2,x,3,y,z,2} or {'df,f,u,2,x,y,z,2} 888 % drops the x-derivative(s) 889 dv:=cddr der; 890 while dv do << 891 if car dv=x then if cdr dv and fixp cadr dv then dv:=cdr dv 892 else 893 else newdv:=cons(car dv,newdv); 894 dv:=cdr dv 895 >>; 896 return if newdv then cons('df,cons(cadr der,reverse newdv)) 897 else cadr der 898end$ 899 900 901symbolic procedure strip_x(ex,ft,x)$ 902begin scalar tm,cex$ 903 % ex is a term 904 % ft is a list of all functions of x which possibly occur in ex 905 return 906 if freeoflist(ex,ft) then 1 else 907 if not pairp ex then ex else 908 if car ex='minus then strip_x(cadr ex,ft,x) else 909 if car ex='df then drop_x_dif(ex,x) else 910 if car ex='expt then if not pairp cadr ex then ex else 911 if caadr ex='df then 912 {'expt,drop_x_dif(cadr ex,x),caddr ex} else 1 % strange 913 else 914 if car ex='times then << 915 ex:=cdr ex; 916 while ex do << 917 cex:=car ex; ex:=cdr ex; 918 if not freeoflist(cex,ft) then 919 if not pairp cex then tm:=cons(cex,tm) else 920 if car cex='df then tm:=cons(drop_x_dif(cex,x),tm) else 921 if car cex='expt then if not pairp cadr cex then tm:=cons(cex,tm) else 922 if caadr cex='df then 923 tm:=cons({'expt,drop_x_dif(cadr cex,x),caddr cex},tm) 924 % else strange - no polynomial in ft 925 >>; 926 if null tm then 1 % strange 927 else 928 if length tm > 1 then reval cons('times,tm) % product of factors 929 else car tm % single factor 930 >> else 1 % strange 931end$ 932 933symbolic procedure sort_partition(pname,p,ft,x)$ 934% The equation is either given by its name pname or its value p 935% if keep_parti=t then the partitioning will be stored 936begin scalar stcp,pcop,parti; 937 % parti will be the list of partial sums 938 if pname then 939 if get(pname,'partitioned) then return get(pname,'partitioned) 940 else <<cp_sq2p_val(pname)$p:=get(pname,'pval)>>$ 941 if (not pairp p) or (car p neq 'plus) then p:=list p 942 else p:=cdr p; 943 if null ft then parti:={{1,1,p}} else 944 while p do << % sort each term into a partial sum 945 stcp:=strip_x(car p,ft,x); % first strip off x_dependent details 946 pcop:=parti; % search for the label in parti 947 while pcop and 948 caar pcop neq stcp do pcop:=cdr pcop; 949 if null pcop then parti:=cons({stcp,1,{car p}},parti) 950 % open a new partial sum 951 else rplaca(pcop,{stcp,add1 cadar pcop, 952 cons(car p,caddar pcop)}); 953 % add the term to an existing partial sum 954 p:=cdr p 955 >>; 956 if pname and keep_parti then put(pname,'partitioned,parti)$ 957 return parti 958end$ % of sort_partition 959 960symbolic procedure int_partition(p,ft,x)$ 961begin scalar parti,ft,pcop; 962 963 % the special case of a quotient 964 if (pairp p) and (car p='quotient) then return 965 if not freeoflist(caddr p,ft) then list p 966 else << 967 pcop:=int_partition(cadr p,ft,x)$ 968 for each h in pcop collect {'quotient,h,caddr p} 969 >>$ 970 parti:=sort_partition(nil,p,ft,x)$ 971 parti:=idx_sort for each h in parti collect cdr h; 972 return for each h in parti collect 973 if car h = 1 then caadr h 974 else cons('plus,cadr h) 975 976end$ 977 978symbolic procedure intpde_(p,ftem,vl,x,potint)$ 979% integration of an polynomial expr. p w.r.t. x 980begin scalar f,ft,l,l1,l2,l3,l4,vl,k,s,a,iflag,flag$ 981 ft:=ftem$ 982 vl:=cons(x,delete(x,vl))$ 983 while ftem and not flag do 984 <<f:=car ftem$ % integrating all terms including car ftem 985 if member(x,fctargs f) then 986 <<if (pairp p) and (car p = 'quotient) then << 987 l1:=lderiv(cadr p,f,vl)$ % numerator 988 if cdr l1 neq 'infinity then << 989 l2:=lderiv(caddr p,f,vl)$ % denomiator 990 if cdr l2 = 'infinity then l1:=l2 991 else << 992 if car l1 and (car l1=car l2) then l1:=(car l1 . 2) % nonlinearity 993 else << 994 l:=lderiv({'plus,car l1,car l2},f,vl)$ % comparison of both 995 if car l2 % changed from car l1 as car l1 gives endless loops with 996 % call intpde((quotient (minus 1) (df u y)),(u),(y x t),y,nil) 997 and (car l = car l2) then l1:=(car l2 . -1) 998 >> 999 >> 1000 >>; 1001 l:=l1; 1002 >> else 1003 l1:=l:=lderiv(p,f,vl)$ 1004 while not (flag or << 1005 iflag:=intlintest(l,x); 1006 if (iflag='noxdrv) or (iflag='nodrv) then << 1007 l2:=start_let_rules()$ 1008 p:=reval p$ 1009 stop_let_rules(l2)$ 1010 l:=lderiv(p,f,vl)$ 1011 iflag:=intlintest(l,x) 1012 >> else 1013 if potint and (iflag='nonlin) and (pairp p) and 1014 (car p='plus) and null l3 and (cdr l neq 'infinity) then << 1015 l3:=t; l4:=l; 1016 l:=car l . 1; 1017 iflag:=intlintest(l,x); 1018 l:=l4 1019 >>; 1020 iflag 1021 >> ) do 1022 <<if (pairp p) and (car p = 'quotient) then 1023 k:=reval {'quotient,coeffn(cadr p,car l,cdr l),caddr p} 1024 else 1025 k:=reval coeffn(p,car l,cdr l)$ 1026 if intcoefftest(car lderiv(k,f,vl),car l,vl) then 1027 <<a:=decderiv(car l,x)$ 1028 %k:=reval list('int,subst('v_a_r_,a,k),'v_a_r_)$ 1029 k:=err_catch_int(subst('v_a_r_,a,k),'v_a_r_)$ 1030 if null k then <<k:=0;flag:='too_slow>>; 1031 if lin_problem then l4:=nil 1032 else l4:=zero_den(k,ft); 1033 if l4 then << 1034 % no constants of integration have been added yet so 1035 % added cases do not involve constants of integration 1036 % for each l2 in l4 do 1037 % to_do_list:=union(list(list('split_into_cases,l2)),to_do_list); 1038 1039 % The above to_do_list:=.. statement can only be issued if 1040 % that expression ex is known not to be an equation or 1041 % equal zero modulo equations. Reason: If the equation 1042 % ex=0 exists then the case ex <> 0 is quickly lead to a 1043 % conradiction and the case ex=0 does not give anything 1044 % new and will execute again the same integration and the 1045 % same case distinction. 1046 1047 flag:='needs_case_split 1048 >> else << 1049 k:=reval subst(a,'v_a_r_,k)$ 1050 s:=cons(k,s)$ 1051 p:=list('difference,p,list('df,k,x))$ %##### This can take long to simplify 1052 p:=err_catch_reval p$ 1053 if null p then flag:='reval_too_slow else 1054 if diffrelp(l1,(l:=lderiv(p,f,vl)),vl) then flag:='neverending 1055 else l1:=l 1056 >> 1057 >> else 1058 flag:='coeffld 1059 >>$ 1060 % iflag='nofct is the so far integrable case 1061 % the non-integrable cases are: flag neq nil, 1062 % (iflag neq nil) and (iflag neq 'nofct), an exception to the 1063 % second case is potint where non-integrable rests are allowed 1064 if (flag neq 'reval_too_slow) and 1065 (flag neq 'neverending) then % because in these cases p is lost, 1066 % it would need to be programmed that 1067 % a trial integration is done. 1068 % no time for doing that now (19.4.2014) 1069 if iflag='nofct then ftem:=smemberl(ftem,p) 1070 else if potint or (fctlength f<length vl) then 1071 <<ftem:=cdr ftem$flag:=nil>> else 1072 flag:=(iflag or flag) 1073 >> else 1074 ftem:=cdr ftem 1075 >>$ 1076 return 1077 if flag then nil else 1078 << 1079 l:=explicitpart(p,ft,x)$ 1080 %l1:=list('INT,l,x)$ 1081 l1:=err_catch_int(l,x)$ 1082 if null l1 then nil 1083 else << 1084 s:=reval cons('plus,cons(l1,s))$ 1085 if lin_problem then l4:=nil 1086 else l4:=zero_den(s,ft); 1087 if l4 then << 1088 % no constants of integration have been added yet so 1089 % added cases do not involve constants of integration 1090 for each l2 in l4 do 1091 % <<write"CCCC"$terpri()$ 1092 to_do_list:=union(list(list('split_into_cases,l2)),to_do_list); 1093 % >>$ 1094 nil 1095 >> else 1096 if freeint_ and (null freeof(s,'INT)) then nil else 1097 if freeabs_ and (null freeof(s,'ABS)) then nil else << 1098 k:=start_let_rules()$ 1099 l2 := reval {'df,l1,x} where !*precise=nil; 1100 if 0 neq (reval {'difference,l,l2} where !*precise=nil) then << 1101 write"REDUCE integrator error:"$terpri()$ 1102 algebraic write "int(",l,",",x,") neq ",l1;terpri()$ 1103 write"Result ignored.";terpri()$ 1104 stop_let_rules(k)$ 1105 nil 1106 >> else << 1107 p:=reval list('difference,p,l2)$ 1108 stop_let_rules(k)$ 1109 if poly_only then if ratexp(s,ft) then list(s,p) 1110 else nil 1111 else list(s,p) 1112 >> 1113 >> 1114 >> 1115 >> 1116end$ % of intpde_ 1117 1118symbolic procedure explicitpart(p,ft,x)$ 1119% part of a sum p which only explicitly depends on a variable x 1120begin scalar l$ 1121 if not member(x,argset smemberl(ft,p)) then l:=p 1122 else if pairp p then 1123 <<if car p='minus then l:=reval list('minus,explicitpart(cadr p,ft,x))$ 1124 if (car p='quotient) and not member(x,argset smemberl(ft,caddr p)) then 1125 l:=reval list('quotient,explicitpart(cadr p,ft,x),caddr p)$ 1126 if car p='plus then << 1127 for each a in cdr p do 1128 if not member(x,argset smemberl(ft,a)) then l:=cons(a,l)$ 1129 if pairp l then l:=reval cons('plus,l)>> >>$ 1130 if not l then l:=0$ 1131 return l$ 1132end$ 1133 1134symbolic procedure intconst(co,vl,x,ifalist)$ 1135% The factors in ifalist must be in the order of integrations done 1136if null ifalist then 0 else 1137begin scalar l,l2,f,coli,cotmp$ 1138 while ifalist do << 1139 cotmp:=coli; 1140 coli:=nil; 1141 while cotmp do << 1142 coli:=cons(list('int,list('times,car ifalist,car cotmp),x),coli); 1143 cotmp:=cdr cotmp 1144 >>; 1145 coli:=cons(1,coli); 1146 ifalist:=cdr ifalist 1147 >>; 1148 1149 while coli do 1150 <<f:=newfct(fname_,delete(x,vl),nfct_)$ 1151 nfct_:=add1 nfct_$ 1152 fnew_:=cons(f,fnew_)$ 1153 flin_:=fctinsert(f,flin_)$ 1154 l:=cons(list('times,f,car coli),l)$ 1155 coli:=cdr coli 1156 >>$ 1157 if length l>1 then l:=cons('plus,l) 1158 else if pairp l then l:=car l 1159 else l:=0$ 1160 if co and co neq 1 then 1161 if pairp co then 1162 <<if car co='times then co:=cdr co 1163 else co:=list co$ 1164 while co do <<if my_freeof(car co,x) then l2:=cons(car co,l2)$ 1165 co:=cdr co>> 1166 >> 1167 else if co neq x then l2:=list co$ 1168 return reval if l2 then cons('times,cons(l,l2)) 1169 else l 1170end$ 1171 1172symbolic procedure intcoefftest(l,l1,vl)$ 1173begin scalar s$ 1174return if not pairp l then t 1175 else if car l='df then 1176 <<s:=decderiv(l1,car vl)$ 1177 if pairp s and pairp cdr s then s:=cddr s 1178 else s:=nil$ 1179 if diffrelp(cons(cddr l,1),cons(s,1),vl) then t 1180 else nil>> 1181else t$ 1182end$ 1183 1184symbolic procedure fctlinear(p,f)$ 1185<<p:=ldiffp(p,f)$ 1186(null car p) and (cdr p=1)>>$ 1187 1188symbolic procedure intlintest(l,x)$ 1189% Test,ob in einem Paar (fuehrende Ableitung.Potenz) 1190% eine x-Ableitung linear auftritt 1191if not car l then 1192 if zerop cdr l then 'nofct 1193 else 'nonlin 1194else % Fkt. tritt auf 1195 if (car l) and (cdr l=1) then % fuehr. Abl. tritt linear auf 1196 if pairp car l then % Abl. der Fkt. tritt auf 1197 if caar l='df then 1198 if member(x,cddar l) then nil 1199 % x-Abl. tritt auf 1200 else if member(x,fctargs cadar l) then 'noxdrv 1201 else 'noxdep 1202 else 'nodrv 1203 else 'nodrv 1204 else 'nonlin$ 1205 1206symbolic procedure decderiv(l,x)$ 1207% in Diff.ausdr. l wird Ordn. d. Abl. nach x um 1 gesenkt 1208begin scalar l1$ 1209return if l then if car l='df then 1210 <<l1:=decderiv1(cddr l,x)$ 1211 if l1 then cons('df,cons(cadr l,l1)) 1212 else cadr l>> 1213 else l 1214 else nil$ 1215end$ 1216 1217symbolic procedure decderiv1(l,x)$ 1218if null l then nil 1219else 1220if car l=x then 1221 if cdr l then 1222 if numberp cadr l then 1223 if cadr l>2 then cons(car l,cons(sub1 cadr l,cddr l)) 1224 else cons(car l,cddr l) 1225 else cdr l 1226 else nil 1227else cons(car l,decderiv1(cdr l,x))$ 1228 1229symbolic procedure integratede(q,ftem,genflag)$ 1230% Integration of a de 1231% result: newde if successfull, nil otherwise 1232begin scalar l,l1,l2,fl,ltdl$ 1233 ftem:=smemberl(ftem,q)$ 1234 ltdl:=length to_do_list$ 1235 1236 again: 1237 if l1:=integrableode(q,ftem) then % looking for an integrable ode 1238 if l1:=integrateode(q,car l1,cadr l1,caddr l1,ftem) then 1239 % trying to integrate it 1240 <<l:=append(cdr l1,l); 1241 q:=prepsq car simplifypdeSQ(simp car l1,ftem,nil,nil,nil)$ 1242 ftem:=smemberl(union(fnew_,ftem),q)$ 1243 fl:=t 1244 >> else 1245 if (ltdl < length to_do_list) and 1246 (caar to_do_list = 'split_into_cases) then 1247 return nil$ % because the ODE involves parameters and the solution 1248 % depends on the value of parameters leading to case distinctions. 1249 % Continuing the integration with other methods would most likely 1250 % lead to the same case distinctions which are to be made next. 1251 1252 if l1:=integratepde(q,ftem,nil,genflag,potint_) then 1253 % trying to integrate a pde 1254 <<q:=car l1$ 1255 for each a in cdr l1 do 1256 <<ftem:=union(fnew_,ftem)$ 1257 if (l2:=integratede(a,ftem,nil)) then l:=append(l2,l) 1258 else l:=cons(a,l) 1259 >>$ 1260 fl:=t$ 1261 if null genflag then l1:=nil$ 1262 ftem:=smemberl(union(fnew_,ftem),q); 1263 goto again 1264 >>$ 1265 1266 if fl then 1267 <<l:=cons(q,l)$ 1268 l:=for each a in l collect reval a$ 1269 l:=for each a in l collect 1270 if pairp a and (car a='quotient) then cadr a 1271 else a>>$ 1272 return l$ 1273end$ 1274 1275symbolic procedure intflagtest(q,fullint)$ 1276if flagp(q,'to_int) then 1277 if fullint then 1278 if (null flagp(q,'to_fullint)) then nil else 1279 if get(q,'starde) then nil else 1280 if (null get(q,'allvarfcts)) 1281 % or (cdr get(q,'allvarfcts)) % if more than one allvar-function 1282 then nil else 1283 begin scalar fl,vl,dl,l,n,mi$ 1284 n:=get(q,'nvars)$ 1285 for each f in get(q,'rational) do % only rational fcts 1286 if fctlength f=n then fl:=cons(f,fl)$ 1287 if null fl then return nil$ 1288 vl:=get(q,'vars)$ 1289 dl:=get(q,'derivs)$ 1290 for each d in dl do 1291 if member(caar d,fl) then 1292 put(caar d,'maxderivs,maxderivs(get(caar d,'maxderivs),cdar d,vl))$ 1293 dl:=fl$ 1294 while vl do 1295 <<mi:=car get(car fl,'maxderivs)$ 1296 l:=list car fl$ 1297 put(car fl,'maxderivs,cdr get(car fl,'maxderivs))$ 1298 for each f in cdr fl do 1299 <<if (n:=car get(f,'maxderivs))=mi then l:=cons(f,l) 1300 else if n<mi then 1301 <<l:=list f$ 1302 mi:=n>>$ 1303 put(f,'maxderivs,cdr get(f,'maxderivs)) 1304 >>$ 1305 dl:=intersection(l,dl)$ 1306 if dl then vl:=cdr vl 1307 else vl:=nil>>$ 1308 for each f in fl do remprop(f,'maxderivs)$ 1309 if fullint and (null dl) then remflag1(q,'to_fullint)$ 1310 return dl 1311 end 1312 else t$ 1313 1314 1315symbolic procedure integrate(q,genintflag,fullint,pdes)$ 1316% integrate pde q; if genintflag is not nil then indirect 1317% integration is allowed 1318% if fullint is not nil then only full integration is allowed which 1319% in addition leads to a substitution 1320% Currently not used if functions occur 1321% Currently not used:Es wird noch nicht ausgenutzt: 1322% 1) functions that occur rationally, 1323% 2) starde 1324% parameter pdes only for drop_pde_from_idties(), drop when pdes_ global 1325% and for mkeqSQlist() for adding inequalities 1326begin scalar l,fli,fnew_old,h,loftolist$ 1327 if fli:=intflagtest(q,fullint) then 1328 <<if fullint then <<fnew_old:=fnew_;fnew_:=nil>>$ 1329 cp_sq2p_val(q)$ 1330 loftolist:=length to_do_list$ 1331 if (l:=integratede(get(q,'pval),get(q,'fcts),genintflag)) then << 1332 if fullint then 1333 while fli and 1334 (not null car (h:=ldiffp(car l,car fli)) or 1335 cdr h neq 1 or << 1336 h:=coeffn(car l,car fli,1); 1337 if domainp h then nil 1338 else << 1339 h:=simplifySQ(cadr h,get(q,'fcts),t,nil,t)$ 1340 if h={(1 . 1)} then t else nil 1341 >> 1342 >> 1343 ) do fli:=cdr fli; 1344 if null fli then << 1345 remflag1(q,'to_fullint); 1346 for each f in fnew_ do drop_fct(f)$ 1347 fnew_:=fnew_old; 1348 l:=nil; 1349 if print_ then << 1350 terpri()$write"Not enough integrations to solve for a function"$ 1351 if null lin_problem then << 1352 write", or,"$terpri()$write"substitution prevented through non-linearity." 1353 >> else write"." 1354 >> 1355 >> else 1356 <<fnew_:=union(fnew_old,fnew_)$ 1357 for each f in fnew_ do << 1358 ftem_:=fctinsert(f,ftem_)$ 1359 flin_:=cons(f,flin_) 1360 >>$ 1361 flin_:=sort_according_to(flin_,ftem_); 1362 fnew_:=nil$ 1363 1364 % Is the old equation to be kept because the new integrated (non-linear) 1365 % equation is sufficient but not necessary? (possible outcome of odeconvert) 1366 h:=cdr l; 1367 while h and (car h neq get(q,'pval)) do h:=cdr h; 1368 if h then << % equation q is to be kept 1369 l:=delete(get(q,'pval),l); 1370 l:=cons(q,mkeqSQlist(nil,nil,l,ftem_,get(q,'vars), 1371 allflags_,t,get(q,'orderings),pdes))$ 1372 if print_ then << 1373 terpri()$ 1374 if l and cdr l and cddr l then << 1375 write"The equation ",q," is kept and an additional sufficient"$terpri()$ 1376 write"integral with conditions ",cdr l," are added." 1377 >> else << 1378 write"The equation ",q," is kept and an additional sufficient"$terpri()$ 1379 write"integral ",cadr l," is added." 1380 >>$ 1381 terpri() 1382 >> 1383 >> else << 1384 flag1(q,'to_eval)$ 1385 updateSQ(q,nil,nil,car l,ftem_,get(q,'vars),t,list(0),nil)$ 1386 drop_pde_from_idties(q,pdes,nil)$ 1387 drop_pde_from_properties(q,pdes)$ 1388 l:=cons(q,mkeqSQlist(nil,nil,cdr l,ftem_,get(q,'vars), 1389 allflags_,t,get(q,'orderings),pdes))$ 1390 put(q,'dec_with,nil); % 23.3.99 added --> cycling? 1391 put(q,'dec_with_rl,nil); % " added --> cycling? 1392 if print_ then << 1393 terpri()$ 1394 if cdr l then 1395 if get(q,'nvars)=get(cadr l,'nvars) then 1396 write "Potential integration of ",q," yields ",l else 1397 write "Partially potential integration of ",q," yields ",l 1398 else write "Integration of ",q$ 1399 terpri() 1400 >> 1401 >>$ 1402 1403 % In order not to re-investigate q for integrability either 1404 % new cases have been added or the integrability is decided: 1405 if loftolist=length to_do_list then << 1406 remflag1(q,'to_fullint)$ 1407 remflag1(q,'to_int) 1408 >> 1409 >> 1410 >> else << 1411 remflag1(q,'to_fullint)$ 1412 remflag1(q,'to_int) 1413 >> 1414 >>$ 1415 return l$ 1416end$ 1417 1418symbolic procedure quick_integrate_one_pde(pdes)$ 1419begin scalar q,p,r,v,nv,nvc,minv,miordr,minofu,minodv,ordr,nofu,nodv$ % ,nvmax$ 1420 % nvmax:=0; 1421 % for each q in ftem_ do if (r:=fctlength q)>nvmax then nvmax:=r; 1422 1423 nv:=no_fnc_of_v()$ % the number of functions for each variable 1424 1425 % find the lowest order derivative wrt. only one variable 1426 miordr:=10000; % the order of the currently best equation 1427 minofu:=10000; % the number of functions depending on the 1428 % variable wrt. which shall be integrated 1429 minodv:=10000; % the number of differentiation variables of 1430 % the so far best equation 1431 while pdes and 1432 (get(car pdes,'length) = 1) do << % only 1 term 1433 q:=get(car pdes,'derivs)$ 1434 if q and % (get(car pdes,'nvars) = nvmax) 1435 cdaar q % any differentiations at all 1436 then << 1437 q:=caar q$ 1438 nodv:=0$ % no of differentiation variables 1439 ordr:=0$ % total order of the derivative 1440 r:=cdr q; 1441 v:=cadr q$ 1442 while r do << 1443 if fixp car r then ordr:=ordr-1+car r 1444 else <<ordr:=add1 ordr;nodv:=add1 nodv>>; 1445 r:=cdr r 1446 >>$ 1447 if nodv>1 then nofu:=10000 % nodv = no of functions depending 1448 else << % on the integration variable 1449 nvc:=nv; 1450 while v neq caar nvc do nvc:=cdr nvc; 1451 nofu:=cdar nvc; 1452 >>$ % no of fncs of v 1453 1454 if nodv=1 then 1455 if ((ordr=1) and (miordr>1)) or 1456 (a_before_b_according_to_c(v,minv,vl_) and (ordr<=miordr)) or 1457 ((v=minv) and 1458 (ordr<miordr) or ((ordr=miordr) and 1459 (nodv<minodv) or ((nodv=minodv) and 1460 (nofu<minofu)))) then 1461 <<p:=car pdes; 1462 minv:=v; 1463 minofu:=nofu; 1464 miordr:=ordr; 1465 minodv:=nodv 1466 >>; 1467 >>$ 1468 pdes:=cdr pdes 1469 >>$ 1470 if p then p:=integrate(p,nil,t,pdes)$ 1471 return p 1472end$ 1473 1474symbolic procedure integrate_one_pde(pdes,genintflag,fullint)$ 1475% trying to integrate one pde 1476begin scalar l,l1,m,p,pdescp,tdlcp$ % ,nvmax,h,f$ 1477 1478 % For non-linear differential equations it can happen, that integrations 1479 % of, for example inhomogeneous term lead to case distinctions on 1480 % constants/functions which have been the result of earlier integrations. 1481 % But if the integration is afterall not completed and therefore not done, 1482 % and therefore no additions should have been added to to_do_list which 1483 % involve only new constants/functions. 1484 tdlcp:=to_do_list$ 1485 1486 % nvmax:=0; 1487 % for each f in ftem_ do if (h:=fctlength f)>nvmax then nvmax:=h; 1488 % at first selecting all eligible de's 1489 m:=-1$ 1490 pdescp:=pdes$ 1491 while pdescp do << 1492 if flagp(car pdescp,'to_int) and not(get(car pdescp,'starde)) and 1493 (null fullint or lin_problem or 1494 get(car pdescp,'linear_) or 1495% not pairp get(car pdescp,'val) or 1496% car get(car pdescp,'val) neq 'times 1497 not pairp get(car pdescp,'fac) 1498 ) then << 1499 l:=cons(car pdescp,l); 1500 if get(car pdescp,'nvars)>m then m:=get(car pdescp,'nvars)$ 1501 >>; 1502 pdescp:=cdr pdescp 1503 >>; 1504 1505 l:=reverse l; 1506 1507 if mem_eff then % try the shortest equation first 1508 while l do 1509 if p:=integrate(car l,genintflag,fullint,pdes) then l:=nil 1510 else l:=cdr l 1511 else 1512 % find an equation to be integrated with as many as possible variables 1513 % if (m=nvmax) or (null fullint) then 1514 while m>=0 do << 1515 l1:=l$ 1516 while l1 do 1517 if (get(car l1,'nvars)=m) and 1518 (p:=integrate(car l1,genintflag,fullint,pdes)) then << 1519 m:=-1$ 1520 l1:=nil 1521 >> else l1:=cdr l1$ 1522 % if fullint then m:=-1 else 1523 m:=sub1 m 1524 >>$ 1525 1526 if null p then << 1527 to_do_list:=cons(1,to_do_list); 1528 L:=to_do_list; 1529 while cdr l neq tdlcp do 1530 if freeoflist(cadr l,ftem_) then rplacd(l,cddr l) 1531 else l:=cdr l; 1532 to_do_list:=cdr to_do_list$ 1533 >>$ 1534return p$ 1535end$ 1536 1537endmodule$ 1538 1539%******************************************************************** 1540module generalized_integration$ 1541%******************************************************************** 1542% Routines for generalized integration of pde's containing unknown 1543% functions 1544% Author: Andreas Brand 1545% December 1991 1546 1547symbolic procedure gintorder(p,vl,x)$ 1548%symbolic procedure gintorder(p,ftem,vl,x)$ 1549% reorder equation p 1550begin scalar l,l1,q,m,b,c,q1,q2$ 1551 if pairp(p) and (car p='quotient) then << 1552 q:=caddr p$ 1553 p:=cadr p$ 1554 l1:=gintorder1(q,x,t)$ 1555% l1:=gintorder1(q,ftem,x,t)$ 1556% if DepOnAllVars(car l1,x,vl) then return nil; 1557 q1:=car l1; 1558 q2:=cadr l1; 1559 >>$ 1560 if pairp(p) and (car p='plus) then p:=cdr p % list of summands 1561 else p:=list p$ 1562 while p do 1563 <<l1:=gintorder1(car p,x,nil)$ 1564% <<l1:=gintorder1(car p,ftem,x,nil)$ 1565 if DepOnAllVars(if q1=1 then car l1 1566 else cons('times, 1567 append(if pairp q1 and car q1='times then cdr q1 1568 else list q1, 1569 if pairp car l1 and caar l1='times then cdar l1 1570 else list car l1)), 1571 x,vl) then l:=p:=nil 1572 else <<l:=termsort(l,l1)$p:=cdr p>> >>$ 1573 if l then 1574 <<l:=for each a in l collect if cddr a then 1575 <<b:=car a$ 1576 c:=cdr reval coeff1(cons('plus,cdr a),x,nil)$ 1577 m:=0$ 1578 while c and (car c=0) do <<c:=cdr c$m:=add1 m>>$ 1579 if m>0 then b:=list('times,list('expt,x,m),b)$ 1580 cons(reval b,c)>> 1581 else 1582 cons(reval car a,cdr reval coeff1(cadr a,x,nil))$ 1583 if q then << 1584 l:=for each a in l collect 1585 cons(car a,for each s in cdr a collect 1586 reval list('quotient,s,q2))$ 1587 l:=for each a in l collect 1588 cons(reval list('quotient,car a,q1),cdr a) 1589 >>$ 1590>>$ 1591return l$ 1592end$ 1593 1594symbolic procedure DepOnAllVars(c,x,vl)$ 1595% tests for occurence of vars from vl in factors of c depending on x 1596begin scalar l$ 1597if pairp c and (car c='times) then c:=cdr c 1598 else c:=list c$ 1599while c and vl do 1600<<if not my_freeof(car c,x) then 1601 for each v in vl do if not my_freeof(car c,v) then l:=cons(v,l)$ 1602 vl:=setdiff(vl,l)$ 1603 c:=cdr c 1604>>$ 1605return (null vl)$ 1606end$ 1607 1608symbolic procedure gintorder1(p,x,mode2)$ 1609%symbolic procedure gintorder1(p,ftem,x,mode2)$ 1610% reorder a term p 1611begin scalar l1,l2,sig$ 1612% mode2 = nil then 1613% l2: list of factors of p not depending on x 1614% or being a power of x 1615% l1: all other factors 1616% mode2 = t then 1617% l2: list of factors of p not depending on x 1618% l1: all other factors 1619 1620if pairp p and (car p='minus) then <<sig:=t$p:=cadr p>>$ 1621if pairp p and (car p='times) then p:=cdr p 1622 else p:=list p$ 1623for each a in p do 1624 <<if my_freeof(a,x) then l2:=cons(a,l2) 1625 % 14 April 2013: dropped 'and freeoflist(a,ftem)' in above if statement 1626 % because then new functions introduced in generalized integrations 1627 % may depend on all variables if all functions in the factors not 1628 % depending on x may depend on all variables apart of x 1629 % 1630 % Example: 0=a(x,y)*b(y,z)-df(c(x,y,z),z) 1631 % If there is 'and freeoflist(a,ftem)' then integration gives 1632 % 0=h(x,y,z) -c(x,y,z)+c_1(x,y), df(h,z)=a*b 1633 % which is useless as the new function has all variables 1634 % and without 'and freeoflist(a,ftem)' integration gives 1635 % 0=a(x,y)*h(y,z)-c(x,y,z)+c_1(x,y), df(h,z)=b 1636 % which is useful as the new function has less variables 1637 % 1638 else if mode2 then l1:=cons(a,l1) 1639 else if a=x then l2:=cons(a,l2) 1640 else if pairp a and (car a='expt) and (cadr a=x) and fixp caddr a 1641 then l2:=cons(a,l2) 1642 else l1:=cons(a,l1)>>$ 1643if pairp l1 then 1644 if cdr l1 then l1:=cons('times,l1) 1645 else l1:=car l1$ 1646if pairp l2 then 1647 if cdr l2 then l2:=cons('times,l2) 1648 else l2:=car l2$ 1649if sig then if l2 then l2:=list('minus,l2) 1650 else l2:=list('minus,1)$ 1651return list(if l1 then l1 else 1,if l2 then l2 else 1)$ 1652end$ 1653 1654symbolic procedure partint(p,ftem,vl,x,genint)$ 1655% genint is the maximal number of terms to be generalized integrated 1656begin scalar f,neg,l1,l2,n,k,l,h$ 1657 if tr_genint then << 1658 terpri()$ 1659 write "generalized integration of the unintegrated rest : "$ 1660 eqprint p 1661 >>$ 1662 l:=gintorder(p,vl,x)$ 1663% l:=gintorder(p,ftem,vl,x)$ 1664 % would too many new equations and functions be necessary? 1665 if pairp(l) and (length(l)>genint) then return nil; 1666 l:=for each s in l collect << 1667 h:=varslist(car s,ftem,vl)$ 1668 if h=nil then << 1669 list('times,x,car s,cons('plus,cdr s)) 1670 >> else << 1671 f:=newfct(fname_,h,nfct_)$ 1672 nfct_:=add1 nfct_$ 1673 fnew_:=cons(f,fnew_)$ 1674 flin_:=fctinsert(f,flin_)$ 1675 neg:=t$ 1676 n:=sub1 length cdr s$ 1677 k:=-1$ 1678 if (pairp car s) and (caar s='df) then 1679 <<h:=reval reval list('difference,cadar s,list('df,f,x,add1 n))$ 1680 if not zerop h then l1:=cons(h,l1)$ 1681 l2:=cddar s>> 1682 else 1683 <<h:=signchange reval reval list('difference,car s, 1684 list('df,f,x,add1 n))$ 1685 if not zerop h then l1:=cons(h,l1)$ 1686 l2:=nil>>$ 1687 reval cons('plus, for each sk on cdr s collect 1688 <<neg:=not neg$ 1689 k:=add1 k$ 1690 reval list('times,if neg then -1 else 1, 1691 append(list('df,f,x,n-k),l2), 1692 tailsum(sk,k,x) ) 1693 >> 1694 ) 1695 >> 1696 >>$ 1697 if l then l:=cons(reval cons('plus,l),l1)$ 1698 if tr_genint then 1699 <<terpri()$ 1700 write "result (without constant or function of integration): "$ 1701 if l then << 1702 eqprint car l$ 1703 write "additional equations : "$ 1704 deprint cdr l 1705 >> else write " nil "$ 1706 >>$ 1707 return l$ 1708end$ 1709 1710symbolic procedure tailsum(sk,k,x)$ 1711begin scalar j$ 1712j:=-1$ 1713return reval cons('plus, 1714for each a in sk collect 1715 <<j:=j+1$ 1716 list('times,a,prod(j+1,j+k),list('expt,x,j)) >> ) 1717end$ 1718 1719symbolic procedure prod(m,n)$ 1720if m>n then 1 1721 else for i:=m:n product i$ 1722 1723endmodule$ 1724 1725 1726%******************************************************************** 1727module intfactor$ 1728%******************************************************************** 1729% Routines for finding integrating factors of PDEs 1730% Author: Thomas Wolf 1731% July 1994 1732 1733% The following without factorization --> faster but less general 1734%symbolic procedure fctrs(p,vl,v)$ 1735%begin scalar fl1,fl2,neg; 1736% 1737%write"p=",p; 1738% 1739% if car p='minus then <<neg:=t;p:=cdr p>>$ 1740% return 1741% if not pairp p then if my_freeof(p,v) and (not freeoflist(p,vl)) then 1742% list(p,1,neg) else 1743% list(1,p,neg) 1744% else if car p='plus then list(1,p,neg) 1745% else 1746% if car p='times then 1747% <<for each el in cdr p do if my_freeof(el,v) and (not freeoflist(p,vl)) then 1748% fl1:=cons(el,fl1) else 1749% fl2:=cons(el,fl2); 1750% if pairp fl1 then fl1:=cons('times,fl1); 1751% if pairp fl2 then fl2:=cons('times,fl2); 1752% if not fl1 then fl1:=1; 1753% if not fl2 then fl2:=1; 1754% list(fl1,fl2,neg) 1755% >> else if my_freeof(p,v) and (not freeoflist(p,vl)) then 1756% list(p,1,neg) else 1757% list(1,p,neg) 1758%end$ % of fctrs 1759% 1760 1761symbolic procedure fctrs(p,indep,v)$ 1762begin scalar fl1,fl2; 1763 p:=cdr reval err_catch_fac p; 1764 for each el in p do 1765 if freeoflist(el,indep) and 1766 ((v=nil) or (not my_freeof(el,v))) then fl1:=cons(el,fl1) 1767 else fl2:=cons(el,fl2); 1768 if null fl1 then fl1:=1; 1769 if null fl2 then fl2:=1; 1770 if pairp fl1 then if length fl1 = 1 then fl1:=car fl1 1771 else fl1:=cons('times,fl1); 1772 if pairp fl2 then if length fl2 = 1 then fl2:=car fl2 1773 else fl2:=cons('times,fl2); 1774 return list(fl1,fl2) 1775end$ % of fctrs 1776 1777 1778symbolic procedure extractfac(p,indep,v)$ 1779% looks for factors of p dependent of v and independent of indep 1780% and returns a list of the numerator factors and a list of the 1781% denominator factors 1782begin scalar nu,de$ 1783 return 1784 if (pairp p) and (car p='quotient) then 1785 <<nu:=fctrs( cadr p,indep,v); 1786 de:=fctrs(caddr p,indep,v); 1787 list( reval if car de neq 1 then list('quotient, car nu, car de) 1788 else car nu, 1789 if cadr de neq 1 then list('quotient, cadr nu, cadr de) 1790 else cadr nu 1791 ) 1792 >> else fctrs(p,indep,v) 1793end$ % of extractfac 1794 1795%---------------------------- 1796 1797symbolic procedure get_kernels(ex)$ 1798% gets the list of all kernels in ex 1799begin scalar res,pri$ 1800 % pri:=t; 1801 ex:=reval ex$ 1802 if pri then <<terpri()$write"ex=",ex>>; 1803 if pairp ex then 1804 if (car ex='quotient) or (car ex='plus) or (car ex='times) then 1805 for each s in cdr ex do res:=union(get_kernels s,res) else 1806 if (car ex='minus) or 1807 ((car ex='expt) and 1808% (numberp caddr ex)) then % not for e.g. (quotient,2,3) 1809 (cadr ex neq 'E) and 1810 (cadr ex neq 'e) and 1811 (not fixp cadr ex) ) then res:=get_kernels cadr ex 1812 else res:=list ex 1813 else if idp ex then res:=list ex$ 1814 if pri then <<terpri()$write"res=",res>>; 1815 return res$ 1816end$ 1817 1818%------------------ 1819 1820symbolic procedure specialsol(p,vl,fl,x,indep,gk)$ 1821% tries a power ansatz for the functions in fl in the kernels 1822% of p to make p to zero 1823% indep is a list of kernels, on which the special solution should 1824% not depend. Is useful, to reduce the search-space, e.g. when an 1825% integrating factor for a linear equation for, say f, is to be 1826% found then f itself can not turn up in the integrating factor fl 1827% gk are kernels which occur in p and possibly extra ones which 1828% e.g. are not longer in p because of factorizing p but which are 1829% likely to play a role, if nil then determined below 1830% x is a variable on which each factor in the special solution has 1831% to depend on. 1832begin 1833 scalar e1,e2,n,nl,h,hh,ai,sublist,eqs,startval,pri,printold,pcopy; 1834% pri:=t; 1835 p:=num p; 1836 pcopy:=p; 1837 if pri then << 1838 terpri()$write"The equation for the integrating factor:"; 1839 terpri()$eqprint p; 1840 >>; 1841 if null gk then gk:=get_kernels(p); 1842 for each e1 in fl do << 1843 h:=nil; %---- h is the power ansatz 1844 if pri then 1845 for each e2 in gk do << 1846 terpri()$write"e2=",e2; 1847 if my_freeof(e2,x) then write" freeof1"; 1848 if not freeoflist(e2,fl) then write" not freeoflist"$ 1849 if not freeoflist(e2,indep) then write" dependent on indep" 1850 >>; 1851 %----- nl is the list of constants to be found 1852 for each e2 in gk do 1853 if (not my_freeof(e2,x)) and % integ. fac. should depend on x 1854 freeoflist(e2,fl) and % the ansatz for the functions to be 1855 % solved should not include these functions 1856 freeoflist(e2,indep) then << 1857 n:=gensym();nl:=cons(n,nl); 1858 h:=cons(list('expt,e2,n),h); 1859 >>; 1860 if h then << 1861 if length h > 1 then h:=cons('times,h) 1862 else h:=car h; 1863 %-- the list of substitutions for the special ansatz 1864 sublist:=cons((e1 . h),sublist); 1865 if pri then <<terpri()$write"Ansatz: ",e1," = ",h>>; 1866 p:= reval num reval subst(h,e1,p); 1867 if pri then <<terpri()$write"p=";eqprint p>> 1868 >> 1869 >>; 1870 if null h then return nil; 1871 %------- An numerical approach to solve for the constants 1872 if nil then << % numerical approach 1873 %--- Substituting all kernels in p by numbers 1874 on rounded; 1875 precision 20; 1876 terpri()$terpri()$write"Before substituting numerical values:"; 1877 eqprint p; 1878 terpri()$terpri()$write"Constants to be calculated: "; 1879 for each n in nl do write n," "; 1880 1881 for each e1 in nl do << 1882 h:=p; 1883 for each e2 in gk do 1884 if freeoflist(e2,fl) then 1885 if pairp e2 and ((car e2 = 'df) or (car e2 = 'int)) then << 1886 n:=list('quotient,1+random 30028,30029); 1887 terpri();write"substitution done: ",e2," = ",n; 1888 h:=subst(list('quotient,1+random 30028,30029),e2,h) 1889 >>; 1890 for each e2 in gk do 1891 if freeoflist(e2,fl) then 1892 if not(pairp e2 and ((car e2 = 'df) or (car e2 = 'int))) then << 1893 n:=list('quotient,1+random 30028,30029); 1894 terpri();write"substitution done: ",e2," = ",n; 1895 h:=subst(list('quotient,1+random 30028,30029),e2,h) 1896 >>; 1897 1898 terpri()$terpri()$write"The equation after all substitutions: "; 1899 terpri()$ 1900 eqprint h; 1901 1902 eqs:=cons(reval h,eqs); 1903 startval:=cons(list('equal,e1,1),startval) 1904 >>; 1905 if length eqs = 1 then eqs:=cdr eqs 1906 else eqs:=cons('list,eqs); 1907 if length startval = 1 then startval:=cdr startval 1908 else startval:=cons('list,startval); 1909 terpri()$write"start rdsolveeval!";terpri()$terpri()$ 1910 h:=rdsolveeval list(eqs,startval); 1911 eqs:=nil; 1912 off rounded; 1913 >>; 1914 1915 %----- An analytical approach to solve for the constants 1916 if null pri then <<printold:=print_;print_:=nil>>; 1917 if p and not zerop p then % uebernommen aus SEPAR 1918 if not (pairp p and (car p='quotient) and % " " " 1919 intersection(argset smemberl(fl,cadr p),vl)) then 1920 p:=separ2(p,fl,list x) else 1921 p:=nil; 1922 if null pri then print_:=printold; 1923 1924 if p then << 1925 % possibly investigating linear dependencies of different caar p 1926 % solve(lasse x-abhaengige und nl-unabhaengige faktoren fallen von 1927 % factorize df(reval list('quotient, caar p1, caar p2),x),nl). 1928 while p do 1929 if freeoflist(cdar p,nl) then <<eqs:=nil;p:=nil>> 1930 % singular system --> no solution 1931 else << 1932 eqs:=cons(cdar p,eqs); 1933 p:=cdr p 1934 >>; 1935 >>; 1936 if pri then <<terpri()$write"eqs1=",eqs>>; 1937 if (null eqs) or (length eqs > maxalgsys_) then return nil 1938 else << 1939 if pri then << 1940 terpri()$write"The algebraic system to solve for ",nl," is:"; 1941 if length eqs > 1 then deprint eqs 1942 else eqprint car eqs 1943 >>; 1944 if length eqs > 1 then eqs:=cons('list,eqs) 1945 else eqs:=car eqs; 1946 1947 if pri then <<terpri()$write"eqs2=",eqs;terpri();write"nl=",nl>>$ 1948 1949 % for catching the error message `singular equations' 1950 hh:=cons('list,nl); 1951 eqs:=<< 1952 ai:=!!arbint; 1953 err_catch_solve(eqs,hh) 1954 >>; 1955 1956 if pri then <<terpri()$write"eqs3a=",eqs," ai=",ai," !!arbint=", 1957 !!arbint;terpri()>>$ 1958 if not freeof(eqs,'arbcomplex) then << 1959 eqs:=reval car eqs; 1960 for h:=(ai+1):!!arbint do 1961 eqs:=subst(0,list('arbcomplex,h),eqs); 1962 if pri then <<terpri()$write"eqs3b=",eqs;terpri()>>$ 1963 eqs:=err_catch_solve(eqs,hh) 1964 >>; 1965 1966 if pri then <<terpri()$write"eqs3c=",eqs;terpri()>>$ 1967 1968 %--- eqs is the list of solutions for the constant exponents of the 1969 %--- integrating factor 1970 1971 if null eqs then return nil; 1972 if length nl=1 then eqs:=list eqs; 1973 if pri then <<write"nl=",nl," eqs4=",eqs;terpri()>>; 1974 1975 for each e1 in eqs do << % each e1 is a list of substitutions 1976 if pri then <<terpri()$write"e2=",e1;terpri()>>$ 1977 if car e1='list then e1:=cdr e1; 1978 if pri then <<terpri()$write"e3=",e1;terpri()>>$ 1979 for each e2 in e1 do << 1980 if pri then algebraic write"solution:",symbolic e2; 1981 sublist:=subst(caddr e2,cadr e2,sublist) 1982 >>; 1983 if pri then <<terpri()$write"The sublist is:",sublist>> 1984 >>; 1985 >>; 1986 if pri then <<terpri()$write"pcopy1=",pcopy;terpri()>>$ 1987 for each e1 in sublist do << 1988 pcopy:=subst(cdr e1,car e1,pcopy); 1989 if pri then <<terpri()$write"e1=",e1;terpri(); 1990 write"pcopy2=",pcopy;terpri()>>$ 1991 >>$ 1992 if pri then <<terpri()$write"pcopy3=",reval pcopy;terpri()>>$ 1993 if pri then <<terpri()$write"pcopy4=",reval reval pcopy;terpri()>>$ 1994 if not zerop reval reval pcopy then return nil else 1995 return for each e1 in sublist collect (car e1 . reval cdr e1) 1996end$ % of specialsol 1997 1998%------------------ 1999 2000symbolic procedure add_factors(gk,p)$ 2001% gk is a list of factors and p anything but a quotient 2002if null p then gk else 2003if ( not pairp p ) or 2004 (( pairp p ) and 2005 (car p neq 'times) ) then 2006union(list if (pairp p) and (car p = 'minus) then cadr p 2007 else p,gk) else 2008<<for each h in cdr p do 2009 gk:=union(list if (pairp h) and (car h = 'minus) then cadr h 2010 else h,gk); 2011 gk 2012>>$ 2013 2014%------------------ 2015 2016symbolic operator findintfac$ 2017symbolic procedure findintfac(pl,ftem,vl,x,doneintvar,intfacdep, 2018 factr,verbse)$ 2019 2020% - pl is a list of equations from which the *-part (inhomogeneous 2021% terms) have been dropped. 2022% - each equation of pl gets an integrating factor h 2023% - doneintvar is a list of variables, on which the integrating factor 2024% should not depend. The chances to find an integrating factor 2025% increase if the inhomogeneous part of pl is dropped and 2026% separately integrated with generalized integration later. 2027% - if factr is not nil then the equation(s) pl is(are) at first 2028% factorized, e.g. if integration(s) have already been done 2029% and there is a chance that the equation can factorize, thereby 2030% simplify and giving a higher chance for integrability. 2031 2032begin 2033 scalar h,newequ,tozero,fl,e1,pri,factr,exfactors,ftr,gk,precise_old, 2034 carftr,zd; 2035 % exfactors is the list of factors extracted at the beginning 2036 % pri:=t; 2037 2038 precise_old:=!*precise$ 2039 !*precise:=nil$ 2040 2041 factr:=t; % whether tozero should be factorized 2042 if pri then <<terpri()$write"START VON FINDINTFAC">>; 2043 %--- Generation of the condition for the integrating factor(s) in fl 2044 for each e1 in pl do << 2045 %--- extracting factors dependend on x and independent of 2046 %--- doneintvar but only if integrations have already been done, 2047 %--- i.e. (doneintvar neq nil) 2048 gk:=union(gk,get_kernels(e1)); 2049 if factr then <<ftr:=extractfac(e1,append(doneintvar,ftem),x); 2050 carftr:=car ftr; 2051 if not evalnumberp carftr then 2052 if (pairp carftr) and 2053 (car carftr='quotient) then 2054 gk:=add_factors(add_factors(gk,cadr carftr), 2055 caddr carftr ) else 2056 gk:=add_factors(gk,carftr); 2057 >> 2058 else ftr:=list(1,nil); 2059 exfactors:=cons(car ftr,exfactors); 2060 if car ftr neq 1 then << 2061 e1:=cadr ftr; 2062 if pri then <<terpri()$write"extracted factor:",eqprint car ftr>>; 2063 >>; 2064 %--- fl is to become the list of integrating factors h 2065 h:=gensym(); 2066 depl!*:=cons(list(h,x),depl!*)$ 2067 depend h,x; 2068 fl:=h . fl; 2069 e1:=intpde(reval list('times,h,e1),ftem,vl,x,t); 2070 if e1 and car e1 then << 2071 newequ:=car e1 . newequ; 2072 tozero:=cadr e1 . tozero; 2073 if pri then << 2074 terpri()$write" the main part of integration:"$ eqprint(car e1); 2075 terpri()$write"car e1=",car e1; 2076 terpri()$write" the remainder of integration:"$ eqprint(cadr e1); 2077 terpri()$write"cadr e1=",cadr e1; 2078 >> 2079 >>; 2080 >>; 2081 if null tozero then return nil; 2082 %-------- newequ is the integral 2083 newequ:=if length pl > 1 then cons('plus,newequ) 2084 else car newequ; 2085 %-------- tozero is the PDE for the integrating factor 2086 tozero:=reval if length pl > 1 then cons('plus,tozero) 2087 else car tozero; 2088 2089 if pairp tozero and (car tozero='quotient) then tozero:=cadr tozero$ 2090 2091 if factr then << 2092 h:=cdr err_catch_fac(tozero)$ 2093 if h then << % else factorization too long 2094 if pri then <<terpri()$write"The factors of tozero:",h>>; 2095 tozero:=nil; 2096 for each e1 in h do 2097 if smemberl(fl,e1) then tozero:=cons(e1,tozero)$ 2098 tozero:= reval if length tozero > 1 then cons('times,tozero) 2099 else car tozero 2100 >> 2101 >>; 2102 2103 if nil and pri then <<write"tozero =";eqprint tozero >>; 2104 h:=nil; 2105 % actually only those f in ftem, in which pl is nonlinear, but also 2106 % then only integrating factors with a leading derivative low enough 2107 h:=specialsol(tozero,vl,fl,x,append(ftem,doneintvar),gk); 2108 % h:=specialsol(tozero,vl,fl,x,doneintvar,gk); 2109 if pri then <<write"h=",h;terpri()>>; 2110 if h then << 2111 for each e1 in h do << % each e1 is one integrating factor determined 2112 if pri then <<terpri()$write"e1=",e1; 2113 terpri()$write"newequ=",newequ;terpri()>>; 2114 newequ:=reval subst(cdr e1,car e1,newequ); 2115 if pri then <<terpri()$write"newequ=",newequ>>; 2116 >> 2117 >> else if pri then write"no integrating factor"; 2118 2119 %--- delete all dependencies of the functions in fl 2120 %--- must come before the following update 2121 for each e1 in fl do << 2122 depl!*:=delete(assoc(e1,depl!*),depl!*)$ 2123 depl!*:=delete(assoc(mkid(e1,'_),depl!*),depl!*)$ 2124 >>; 2125 2126 %--- update intfacdep 2127 for each e1 in vl do 2128 if (e1 neq x) and my_freeof(intfacdep,e1) and 2129 ((not my_freeof(h,e1)) or (not my_freeof(exfactors,e1))) 2130 then intfacdep:=cons(e1,intfacdep); 2131 2132 %--- returns nil if no integrating factor else a list of the 2133 %--- factors and the integral 2134 if h and print_ and verbse then << 2135 terpri()$write"The integrated equation: "; 2136 eqprint newequ; 2137 terpri()$ 2138 if length pl = 1 then write"An integrating factor has been found:" 2139 else write"Integrating factors have been found: "$ 2140 >>; 2141 !*precise:=precise_old$ 2142 2143 if (null h) or (zerop newequ) then return nil$ 2144 2145 % test for zero denominators 2146 zd:=zero_den(h,ftem); 2147 if zd then return << 2148 % no constants of integration have been added yet so 2149 % added cases do not involve constants of integration 2150 for each e1 in zd do 2151 % <<write"DDDD"$terpri()$ 2152 to_do_list:=union({list('split_into_cases,e1)},to_do_list); 2153 % >>$ 2154 nil 2155 >>$ 2156 2157 return 2158 list(newequ, 2159 for each e1 in h collect << 2160 ftr:=car exfactors; 2161 exfactors:=cdr exfactors; 2162 gk:=if ftr=1 then cdr e1 2163 else reval list('quotient,cdr e1,ftr); 2164 if print_ and verbse then mathprint gk; 2165 gk 2166 >>, 2167 intfacdep) 2168end$ % findintfac 2169 2170endmodule$ 2171 2172 2173%******************************************************************** 2174module odeintegration$ 2175%******************************************************************** 2176% Routines for integration of ode's containing unnown functions 2177% Author: Thomas Wolf 2178% August 1991 2179 2180symbolic procedure integrateode(de,fold,xnew,ordr,ftem)$ 2181begin scalar newde,newnewde,l,h,newcond$ 2182 h:= % the integrated equation 2183 if not xnew then << % Integr. einer alg. Gl. fuer eine Abl. 2184 % The following call of solveeval is taken out as it might 2185 % potentially loose cases for non-linear problems. 2186 % It definitely drops constant ftem-dependent factors which 2187 % correspond to a loss of cases which could be fixed through 2188 % testing : if pairp de and car de = 'times then newde:=de else ... 2189 2190 %newde:=cadr solveeval list(de,fold)$ 2191 %if not freeof(newde,'ROOT_OF) then nil 2192 % else << 2193 % newde:=reval list('plus,cadr newde,list('minus,caddr newde))$ 2194 % if (l:=integratepde(newde,ftem,nil,genint_,nil)) then 2195 2196 if (l:=integratepde(de,ftem,nil,genint_,nil)) then 2197 <<newcond:=append(newcond,cdr l);car l>> 2198 %genflag=t,potflag=nil 2199 else nil 2200 % >> 2201 >> else % eine ode fuer ein f? 2202 if not pairp fold then % i.e. not df(...,...), i.e. fold=f 2203 if (l:=odeconvert(de,fold,xnew,ordr,ftem)) then 2204 <<newcond:=append(newcond,cdr l);car l>> else nil 2205 % --> ode fuer eine Abl. von f 2206 else << 2207 newde:=if (l:=odeconvert(de,fold,xnew,ordr,ftem)) then 2208 <<newcond:=append(newcond,cdr l);car l>> else nil$ 2209 if not newde then nil 2210 else << 2211 % Similarly to above for safety reasons and currently the 2212 % inability to handle more than one solution the following 2213 % solveeval is commented out 2214 2215 % newnewde:=cadr solveeval list(newde,fold)$ 2216 % newnewde:=reval list('plus,cadr newnewde,list('minus, 2217 % caddr newnewde))$ 2218 ftem:=union(fnew_,ftem)$ 2219 % newnewde:=integratede(newnewde,ftem,nil)$ 2220 newnewde:=integratede(newde,ftem,nil)$ 2221 if newnewde then <<newcond:=append(newcond,cdr newnewde); 2222 car newnewde>> 2223 else newde 2224 >> 2225 >>; 2226 2227 return if not h then nil 2228 else cons(h,newcond) 2229 2230end$ % of integrateode 2231 2232symbolic procedure odecheck(ex,fint,ftem)$ 2233% assumes an reval'ed expression ex 2234% Goes wrong if car ex is a list of expressions! 2235begin scalar a,b,op,ex1$ 2236 %***** ex is a ftem-function ***** 2237 if ex=fint then % list(ex,0,0,..) 2238 <<a:=list ex$ 2239 ex:=fctargs ex$ 2240 while ex do 2241 <<a:=append(list(0,0),a)$ 2242 ex:=cdr ex>>$ 2243 % not checked if it is a function of an expression of x 2244 op:=reverse a>> 2245 else if pairp ex then 2246 %***** car ex is 'df ***** 2247 if (car ex)='df then 2248 <<a:=odecheck(cadr ex,fint,ftem)$ 2249 if not pairp a then op:=a 2250 else % a is list(fctname,0,0,..,0,0) 2251 <<op:=list(car a)$ 2252 a:=fctargs car a$ % a is list(variables), not checked 2253 ex:=cddr ex$ % ex is list(derivatives) 2254 while a do 2255 <<ex1:=ex$ 2256 while ex1 and ((car a) neq (car ex1)) do ex1:=cdr ex1$ 2257 if null ex1 then op:=cons(0,cons(0,op)) 2258 else 2259 <<if not cdr ex1 then b:=1 % b is number of deriv. of that var. 2260 else 2261 <<b:=cadr ex1$ 2262 if not numberp b then b:=1>>$ 2263 op:=cons(b,cons(b,op))>>$ 2264 a:=cdr a>>$ 2265 op:=reverse op>> >> 2266 else 2267 %***** car ex is a standard or other function ***** 2268 <<a:=car ex$ % for linearity check 2269 ex:=cdr ex$ 2270 if a='int then ex:=list reval car ex$ 2271 while (op neq '!_abb) and ex do 2272 <<b:=odecheck(car ex,fint,ftem)$ 2273 if b then % function found 2274 if b eq '!_abb then op:='!_abb % occures properly 2275 else op:=odetest(op,b)$ 2276 ex:=cdr ex>> >>$ 2277 return op 2278end$ 2279 2280symbolic procedure integrableode(p,ftem)$ 2281if % length get(p,'derivs) ? 2282 delength p 2283 >(if odesolve_ then odesolve_ else 0) then 2284 (if cont_ then 2285 if yesp("expression to be integrated ? ") then 2286 integrableode1(p,ftem)) 2287else integrableode1(p,ftem)$ 2288 2289symbolic procedure integrableode1(p,ftem)$ 2290begin scalar a,b,u,vl,le,uvar, 2291 fint,fivar,% the function to be integrated and its variables 2292 fold, % the new function of the ode 2293 xnew, % the independ. variable of the ode 2294 ordr1, % order of the ode 2295 ordr2, % order of the derivative for which it is an ode 2296 intlist$ % list of ode's 2297 ftem:=smemberl(ftem,p)$ 2298 vl:=argset ftem$ 2299% p muss genau eine Funktion aus ftem von allen Variablen enthalten. 2300% Die Integrationsvariable darf nicht Argument anderer in p enthaltener 2301% ftem-Funktionen sein. 2302 a:=ftem$ 2303 b:=nil$ 2304 le:=length vl$ 2305 while a and vl do 2306 <<u:=car a$ 2307 uvar:=fctargs u$ 2308 if (length uvar) = le then 2309 if b then 2310 <<vl:=nil$a:=list(nil)>> 2311 else 2312 <<b:=t$ 2313 fint:=u$ 2314 fivar:=uvar>> 2315 else vl:=setdiff(vl,uvar)$ 2316 a:=cdr a>>$ 2317 2318 if not b then vl:=nil$ 2319 le:=length p$ 2320 if ((1<le) and vl) then << 2321 a:=odecheck(p,fint,ftem)$ 2322 if not atom a then << % The equation is an ode 2323 ordr1:=0$ 2324 ordr2:=0$ 2325 xnew:=nil$ 2326 a:=cdr a$ 2327 b:=fivar$ 2328 while b do 2329 <<if (car a) neq 0 then 2330 <<fold:=cons(car b,fold)$ 2331 if (car a) > 1 then fold:=cons(car a,fold)>>$ 2332 ordr2:=ordr2+car a$ 2333 if (car a) neq (cadr a) then 2334 <<xnew:=car b$ 2335 if not member(xnew,vl) then 2336 <<b:=list(nil)$vl:=nil>>$ 2337 ordr1:=(cadr a) - (car a)>>$ 2338 b:=cdr b$ 2339 a:=cddr a>>$ 2340 fold:=reverse fold$ 2341 %fold is the list of diff.variables + number of diff. 2342 if fold then fold:=cons('df,cons(fint,fold)) 2343 else fold:=fint$ 2344 if vl and ((ordr1 neq 0) or (ordr2 neq 0)) then 2345 intlist:=list(fold,xnew,ordr1,ordr2) 2346 >> % of variable found 2347 >>$ % of if 2348 return intlist 2349end$ % of integrableode1 2350 2351symbolic procedure odetest(op,b)$ 2352if not op then b 2353else % op=nil --> first function found 2354 if (car op) neq (car b) then '!_abb else % f occurs in differ. fct.s 2355begin scalar dif,a$ 2356 dif:=nil$ % dif=t --> different derivatives 2357 a:=list(car op)$ % in one variable already found 2358 op:=cdr op$ 2359 b:=cdr b$ 2360 while op do 2361 <<a:=cons(max(cadr op,cadr b),cons(min(car op,car b),a))$ 2362 if (car a) neq ( cadr a) then 2363 if dif then 2364 <<a:='!_abb$ 2365 op:=list(1,1)>> 2366 else dif:=t$ 2367 op:=cddr op$ 2368 b:=cddr b>>$ 2369 if not atom a then a:=reverse a$ 2370 return a % i.e. '!_abb or (fctname,min x1-der.,max x1-der.,...) 2371end$ 2372 2373symbolic procedure odeconvert(de,ford,xnew,ordr,ftem)$ 2374begin scalar j,h,h1,h2,t2,ford_,newco,oldde,newde,newvl,null_,ruli,ld$ 2375% trig1,trig2,trig3,trig4,trig5,trig6,zd$ 2376 ford_:=gensym()$ 2377 depl!*:=delete(assoc(ford_,depl!*),depl!*)$ 2378 depend1(ford_,xnew,t)$ 2379 oldde:=reval subst(ford_,reval ford,de)$ 2380 if pairp ford then % i.e. (car ford) eq 'df then 2381 << for j:=1:ordr do 2382 oldde:= reval subst( reval list('df,ford_,xnew,j), 2383 reval list('df,ford,xnew,j), oldde)>>$ 2384 algebraic !!arbconst:=0$ 2385 %newde:=algebraic(first odesolve(symbolic oldde,symbolic ford_,symbolic xnew))$ 2386 %newde:=algebraic(first lisp err_catch_odesolve(oldde,ford_,xnew))$ 2387 2388 newde:=err_catch_odesolve(oldde,ford_,xnew)$ 2389 if newde and (car newde='list) and cdr newde and cddr newde then return << 2390 if print_ then <<terpri()$ 2391 write "The ode has more than one solution."$ 2392 algebraic write "Equation is: ",algebraic symbolic oldde$ 2393 algebraic write "Solution is: ",algebraic symbolic newde 2394 >>; 2395 nil 2396 >> else 2397 newde:=cadr newde$ 2398 2399 if null newde then return nil; 2400 2401 % Check that the coefficient of the highest power of the highest derivative 2402 % is not free of ftem functions/constants as there can be cases when 2403 % ford_ can have an arbitrary value because all terms with ford_ or a 2404 % derivative of ford_ vanishes and the solution of the ODE does not 2405 % involve quotients or logarithms so special cases would not be found to 2406 % be investigated. 2407 ld:=reval {'df,ford_,xnew,ordr}$ 2408 h:=algebraic(coeff(lisp oldde,lisp ld)); 2409 while cdr h do h:=cdr h$ 2410 h:=simp car h; 2411 h:=simplifySQ(h,ftem,t,nil,t)$ 2412 if null cdr h then h:=car h 2413 else << 2414 j:=car h; h:=cdr h$ 2415 while h do << 2416 j:=multsq(car h,j)$ h:=cdr h 2417 >>$ 2418 h:=j 2419 >>$ 2420 2421 if not freeoflist(h,ftem) and not member(h,ineq_) then return << 2422 % This case can be added because h does not contain new 2423 % constants of integration. 2424 to_do_list:=cons(list('split_into_cases,h),to_do_list); 2425 nil 2426 >>$ 2427 2428 % Check the case that newde has the form 2429 % (equal ford_ (equal .. ..)) 2430 % where .. includes ford_ 2431 % which is an error of ODESOLVE 2432 if pairp newde and 2433 car newde = 'equal and 2434 2435 pairp cdr newde and 2436 cadr newde = ford_ and 2437 2438 pairp cddr newde and 2439 pairp caddr newde and 2440 2441 caaddr newde = 'equal then 2442 2443 if freeof(caddr newde,ford_) then newde:=nil 2444 else << 2445 h:=solveeval {{'list,{'difference,cadr caddr newde,caddr caddr newde}}, 2446 {'list,ford_}}; 2447 if (not freeof(h,'i)) or (not freeof(h,'!:gi!:)) then algebraic(on complex)$ 2448 if pairp h and (car h = 'list) and (pairp cadr h) and 2449 (caadr h = 'equal) and (cadadr h = ford_) then newde:=cadr h 2450 >>$ 2451 2452 % if newde is a parametric solution, i.e. has the form 2453 % (list (equal ford_ ...) (equal xnew ...) arbparam(1)) 2454 if pairp newde and 2455 car newde = 'list and 2456 2457 pairp cdr newde and 2458 pairp cadr newde and 2459 caadr newde = 'equal and 2460 cadadr newde = ford_ and 2461 2462 pairp cddr newde and 2463 pairp caddr newde and 2464 caaddr newde = 'equal and 2465 car cdaddr newde = xnew and 2466 2467 pairp cdddr newde and null 2468 cddddr newde then << 2469 2470 h:=solveeval {{'list,{'difference,cadr caddr newde,caddr caddr newde}}, 2471 {'list,cadddr newde}}; 2472 2473 if (not freeof(h,'i)) or (not freeof(h,'!:gi!:)) then algebraic(on complex)$ 2474 if pairp h and (car h = 'list) and (pairp cadr h) and 2475 (caadr h = 'equal) and (cadadr h = cadddr newde) then 2476 newde:={'equal,ford_,subst(car cddadr h,cadadr h,car cddadr newde)} 2477 >>; 2478 2479 % if the solution is a list of 2 solutions, i.e. a case distinction 2480 % would be necessary then currently these solutions are only printed 2481 % and can not be used. 2482 if pairp newde and 2483 car newde = 'list and 2484 2485 pairp cdr newde and 2486 pairp cadr newde and 2487 caadr newde = 'equal and 2488% cadadr newde = ford_ and 2489 2490 pairp cddr newde and 2491 pairp caddr newde and 2492 caaddr newde = 'equal and 2493% car cdaddr newde = xnew 2494 car cdaddr newde = cadadr newde 2495 2496 then return << 2497 if print_ then <<terpri()$ 2498 write "The ode has more than one solution."$ 2499 algebraic write "Equation is: ",algebraic symbolic oldde$ 2500 algebraic write "Solution is: ",algebraic symbolic newde 2501 >>; 2502 nil 2503 >>$ 2504 2505 ruli:= start_let_rules()$ 2506 2507 if newde then 2508 if !*rational then << off rational; newde:=reval newde; on rational>> 2509 else newde:=reval newde; 2510 2511 % It the solution has the form {'equal,expression,0} then solve the 2512 % numerator of expression for ford_ . 2513 2514 if newde and 2515 (car newde='equal) and 2516 (caddr newde=0) then newde:={'equal,reval num cadr newde,0}; 2517 2518 % Update of flin_ 2519 if flin_ and 2520 newde and 2521 (car newde='equal) then 2522 2523 if not freeoflist(cadr newde,flin_) then << 2524 % The function that is solved for is a non-linearly occuring function 2525 % and therefore all occuring functions are now deleted from flin_. 2526 h1:=flin_; 2527 for each h in h1 do 2528 if not freeof(newde,h) then flin_:=delete(h,flin_) 2529 >> else << 2530 2531 % Are there flin_ constants/functions which do not occure linearly anymore? 2532 j:=nil; 2533 h1:=flin_; 2534 for each h in h1 do 2535 if not freeof(newde,h) then 2536 if null lin_check(caddr newde,{h}) then flin_:=delete(h,flin_) 2537 else j:=cons(h,j); 2538 if j and cdr j then 2539 if null lin_check(caddr newde,j) then 2540 for each h in j do flin_:=delete(h,flin_)$ 2541 >>$ 2542 2543 j:=zero_den(subst(ford,ford_,newde),cons(ford_,ftem)); %,argset ftem 2544 if null j and not freeoflist(caddr newde,ftem_) then 2545 j:=ProportionalityConditions(caddr newde, 2546 for h1:= 1 : (algebraic !!arbconst) 2547 collect list('arbconst,h1), xnew)$ 2548 if j then return 2549 if t then nil 2550 else << 2551 2552 % The following makes only sense if the special case to be investigated 2553 % involves constants of integration. But then these case distinctions 2554 % can not be set up because the the new constants of integration are not 2555 % known after leaving this integration module. 2556 % possible solution: either adding fnew_ to ftem even if no integrated 2557 % equation is returned, or returning a product of the integrated 2558 % equation and the expressions in j. 2559 2560 newvl:=delete(xnew,if (pairp ford and (car ford='df)) 2561 then fctargs cadr ford 2562 else fctargs ford)$ 2563 2564 for h1:=1 : algebraic !!arbconst do 2565 if not freeof(j,list('arbconst,h1)) then << 2566 newco:=newfct(fname_,newvl,nfct_)$ 2567 nfct_:=add1 nfct_$ 2568 fnew_:=cons(newco,fnew_)$ 2569 j:=for each h in j collect subst(newco,list('arbconst,h1),h) 2570 >>$ 2571 2572 while j do << 2573 if freeoflist(car j,fnew_) and not member(car j,cdr j) then 2574 to_do_list:=cons(list('split_into_cases,car j),to_do_list); 2575 j:=cdr j 2576 >>$ 2577 nil 2578 >>$ 2579 2580 % if safeint_ and zero_den(newde,ftem,argset ftem) then newde:=nil; 2581 if freeint_ and null freeof(newde,'int) then newde:=nil; 2582 if freeabs_ and null freeof(newde,'abs) then newde:=nil; 2583 if newde and (cadr newde neq oldde) then << % solution found 2584 % Test der Loesung klappt nur, wenn Loesung explizit gegeben 2585 if cadr newde neq ford_ then << 2586 if print_ then <<terpri()$ 2587 write "Solution of the ode is not explicitly given."$ 2588 algebraic write "Equation is: ",algebraic symbolic oldde$ 2589 algebraic write "Solution is: ",algebraic symbolic newde 2590 >>; 2591 if poly_only then % The solution must be rational in the 2592 % function and constants of integration 2593 if not rationalp(newde,ford_) then newde:=nil else << 2594 j:=1; 2595 while (j leq ordr) and 2596 rationalp(subst(ford_,list('arbconst,j),newde),ford_) do j:=j+1; 2597 if j leq ordr then newde:=nil 2598 >>; 2599 if pairp newde and (car newde = 'equal) then 2600 if (pairp cadr newde) and 2601 (caadr newde = 'quotient) and 2602 (zerop caddr newde) then newde:={'equal,cadadr newde,0} 2603 else 2604 if (pairp caddr newde) and 2605 (caaddr newde = 'quotient) and 2606 (zerop cadr newde) then newde:={'equal,0,cadr caddr newde} 2607 else << 2608 % If it is linear in one arbcons(i) and the second terms in the 2609 % expression is a sqrt or tan or arctan,... then delete this function. 2610 % Example: arbconst(1) + sqrt(g0785**2 + u1**2 )=0 2611 h1:=reval {'difference,cadr newde,caddr newde}; 2612 j:=1; 2613 while (j leq ordr) and 2614 (freeof(h1,list('arbconst,j)) or 2615 (not lin_check(h1,list list('arbconst,j)))) do j:=j+1; 2616 if j leq ordr then << 2617 h:=coeffn(h1,list('arbconst,j),1)$ 2618 t2:=reval {'quotient,{'difference,{'times,list('arbconst,j),h},h1},h}; 2619 h2:={'equal,{'plus,list('arbconst,j), 2620 simplifiziere(t2,cons(ford_,cons(xnew,ftem)),nil)},0}$ 2621 if h2 neq newde then << 2622 newde:=h2$ 2623 if print_ then 2624 algebraic (write "The solution is modified to: ",symbolic newde) 2625 >> 2626 >> 2627 >> 2628 >> else << 2629 null_:=simp!* reval subst(caddr newde, ford_, oldde)$ 2630 if not sqzerop null_ then << 2631 if print_ then << 2632 write "odesolve solves : "$ 2633 deprint list oldde$ 2634 write "to"$ 2635 eqprint newde$ 2636 Write "which inserted in the equation yields "$ 2637 eqprint {'!*sq,null_,t} 2638 >> 2639 >> 2640 >> 2641 >>$ 2642 2643 if newde then 2644 <<newde:=list('plus,cadr newde,list('minus,caddr newde))$ 2645 if zerop reval list('plus,newde,list('minus,oldde)) then newde:=nil$ 2646 % not anymore needed after writing to to_do_list 2647 % if newde and (zd neq nil) then 2648 % newde:=cons('times,append(zd,list newde))$ 2649 >>$ 2650 2651 depl!*:=delete(assoc(ford_,depl!*),depl!*)$ 2652 stop_let_rules(ruli)$ 2653 return 2654 if null newde then nil 2655 else 2656 <<newde:=subst(ford,ford_,newde)$ 2657 newvl:=delete(xnew,if (pairp ford and (car ford='df)) 2658 then fctargs cadr ford 2659 else fctargs ford)$ 2660% if pairp ford then newvl:=delete(xnew,cdr assoc(cadr ford,depl!*)) 2661% else newvl:=delete(xnew,cdr assoc(ford,depl!*))$ 2662 2663 for j:=1 : algebraic !!arbconst do 2664 if not freeof(newde,list('arbconst,j)) then << 2665 newco:=newfct(fname_,newvl,nfct_)$ 2666 nfct_:=add1 nfct_$ 2667 fnew_:=cons(newco,fnew_)$ 2668 newde:=subst(newco,list('arbconst,j),newde)$ 2669 if lin_check(newde,{newco}) then flin_:=fctinsert(newco,flin_)$ 2670 >>$ 2671 {newde} 2672 >> 2673end$ 2674 2675endmodule$ 2676 2677%******************************************************************** 2678module jetdifferentiation$ 2679%******************************************************************** 2680% Routines supporting liepde and conlaw 2681% Author: Thomas Wolf 2682% 1998 2683 2684%------------- 2685 2686symbolic procedure comparedif1(u1l,u2l)$ 2687% checks whether u2l has more or at least equally many 1's, 2's, ... 2688% contained as u1l. 2689% returns a list of 1's, 2's, ... which are in excess in u2l 2690% compared with u1l. The returned value is 0 if both are identical 2691begin 2692 scalar ul; 2693 if u2l=nil then if u1l neq nil then return nil 2694 else return 0 2695 else if u1l=nil then return u2l 2696 else % both are non-nil 2697 if car u1l < car u2l then return nil else 2698 if car u1l = car u2l then return comparedif1(cdr u1l,cdr u2l) else << 2699 ul:=comparedif1(u1l,cdr u2l); 2700 return if not ul then nil else 2701 if zerop ul then list car u2l else 2702 cons(car u2l,ul) 2703 >> 2704end$ % of comparedif1 2705 2706%------------- 2707 2708symbolic procedure comparedif2(u1,u1list,du2)$ 2709% checks whether du2 is a derivative of u1 differentiated 2710% wrt. u1list 2711begin 2712 scalar u2l; 2713 u2l:=combidif(du2)$ % u2l=(u2, 1, 1, ..) 2714 if car u2l neq u1 then return nil else 2715 return comparedif1(u1list, cdr u2l) 2716end$ % of comparedif2 2717 2718%------------- 2719 2720symbolic procedure listdifdif1(du1,deplist)$ 2721% lists all elements of deplist which are *not* derivatives 2722% of du1 2723begin 2724 scalar u1,u1list,res,h$ 2725 h:=combidif(du1); 2726 u1:=car h; 2727 u1list:=cdr h; 2728 for each h in deplist do 2729 if not comparedif2(u1,u1list,h) then res:=cons(h,res); 2730 return res 2731end$ % of listdifdif1 2732 2733%------------- 2734 2735symbolic procedure diffdeg(p,v)$ 2736% liefert Ordnung der Ableitung von p nach v$ 2737% p Liste Varible+Ordnung der Ableitung, v Variable (Atom) 2738if null p then 0 % alle Variable bearbeitet ? 2739else if car p=v then % v naechste Variable ? 2740 if cdr p then 2741 if numberp(cadr p) then cadr p % folgt eine Zahl ? 2742 else 1 2743 else 1 2744 else diffdeg(cdr p,v)$ % weitersuchen 2745 2746%------------- 2747 2748symbolic procedure ldiff1(l,v)$ 2749% liefert Liste der Ordnungen der Ableitungen nach den Variablen aus v 2750% l Liste (Variable + Ordnung)$ v Liste der Variablen 2751if null v then nil % alle Variable abgearbeitet ? 2752else cons(diffdeg(l,car v),ldiff1(l,cdr v))$ 2753 % Ordnung der Ableitung nach 2754 % erster Variable anhaengen 2755 2756%------------- 2757 2758symbolic procedure ldifftot(p,f)$ 2759% leading derivative total degree ordering 2760% liefert Liste der Variablen + Ordnungen mit Potenz 2761% p Ausdruck in LISP - Notation, f Funktion 2762ldifftot1(p,f,fctargs f)$ 2763 2764%------------- 2765 2766symbolic procedure ldifftot1(p,f,vl)$ 2767% liefert Liste der Variablen + Ordnungen mit Potenz 2768% p Ausdruck in LISP - Notation, f Funktion, lv Variablenliste 2769begin scalar a$ 2770 a:=cons(nil,0)$ 2771 if not atom p then 2772% if member(car p,list('expt,'plus,'minus,'times, 2773% 'quotient,'df,'equal)) then 2774 if member(car p,REDUCEFUNCTIONS_) then 2775 % erlaubte Funktionen 2776 <<if (car p='plus) or (car p='times) or 2777 (car p='quotient) or (car p='equal) then 2778 <<p:=cdr p$ 2779 while p do 2780 <<a:=diffreltot(ldifftot1(car p,f,vl),a,vl)$ 2781 p:=cdr p 2782 >> 2783 >> else 2784 if car p='minus then a:=ldifftot1(cadr p,f,vl) else 2785 if car p='expt then % Exponent 2786% if numberp caddr p then 2787% <<a:=ldifftot1(cadr p,f,vl)$ 2788% a:=cons(car a,times(caddr p,cdr a)) 2789% >> else a:=cons(nil,0) 2790 <<a:=ldifftot1(cadr p,f,vl)$ 2791 if (numberp caddr p) and 2792 (numberp cdr a) then a:=cons(car a,times(caddr p,cdr a)) 2793 else a:=cons(car a,10000) 2794 >> 2795 % Potenz aus Basis wird mit 2796 % Potenz multipliziert 2797 else 2798 if car p='df then % Ableitung 2799 if cadr p=f then a:=cons(cddr p,1) 2800 % f wird differenziert? 2801 else a:=cons(nil,0) 2802 else % any other non-linear function 2803 <<p:=cdr p$ 2804 while p do 2805 <<a:=diffreltot(ldifftot1(car p,f,vl),a,vl)$ 2806 p:=cdr p 2807 >>; 2808 a:=cons(car a,10000) 2809 >> 2810 >> else % sonst Konstante bzgl. f 2811 2812 if p=f then a:=cons(nil,1) % Funktion selbst 2813 else a:=cons(nil,0) % alle uebrigen Fkt. werden 2814 else if p=f then a:=cons(nil,1)$ % wie Konstante behandelt 2815 return a 2816end$ 2817 2818%------------- 2819 2820symbolic procedure diffreltot(p,q,v)$ 2821% liefert komplizierteren Differentialausdruck$ 2822if diffreltotp(p,q,v) then q 2823 else p$ 2824 2825%------------- 2826 2827symbolic procedure diffreltotp(p,q,v)$ 2828% liefert t, falls p einfacherer Differentialausdruck, sonst nil 2829% p, q Paare (liste.power), v Liste der Variablen 2830% liste Liste aus Var. und Ordn. der Ableit. in Diff.ausdr., 2831% power Potenz des Differentialausdrucks 2832begin scalar n,m$ 2833m:=eval(cons('plus,ldiff1(car p,v)))$ 2834n:=eval(cons('plus,ldiff1(car q,v)))$ 2835return 2836 if m<n then t 2837 else if n<m then nil 2838 else diffrelp(p,q,v)$ 2839end$ 2840 2841%------------- 2842 2843algebraic procedure subdif1(xlist,ylist,ordr)$ 2844% A list of lists of derivatives of one order for all functions 2845begin 2846 scalar allsub,revx,i,el,oldsub,newsub; 2847 revx:=sqreverse xlist; 2848 allsub:={}; 2849 oldsub:= for each y in ylist collect y=y; 2850 for i:=1:ordr do % i is the order of next substitutions 2851 <<oldsub:=for each el in oldsub join nextdy(revx,xlist,el); 2852 allsub:=sqcons(oldsub,allsub) 2853 >>; 2854 return allsub 2855end$ 2856 2857%------------- 2858 2859algebraic procedure nextdy(revx,xlist,dy)$ 2860% generates all first order derivatives of lhs dy 2861% revx = reverse xlist; xlist is the list of variables; 2862% dy the old derivative 2863begin 2864 scalar x,n,ldy,rdy,ldyx,sublist; 2865 x:=sqfirst revx; revx:=sqrest revx; 2866 sublist:={}; 2867 ldy:=lhs dy; 2868 rdy:=rhs dy; 2869 2870 while lisp(not member(prepsq simp!* algebraic x, 2871 prepsq simp!* algebraic ldy)) 2872 and (revx neq {}) do 2873 <<x:=sqfirst revx; revx:=sqrest revx>>; 2874 2875 n:=length xlist; 2876 if revx neq {} then % dy is not the function itself 2877 while sqfirst xlist neq x do xlist:=sqrest xlist; 2878 xlist:=reverse xlist; 2879 2880 % New higher derivatives 2881 while xlist neq {} do 2882 <<x:=sqfirst xlist; 2883 ldyx:=df(ldy,x); 2884 sublist:=cons((lisp reval algebraic ldyx)= 2885 mkid(mkid(rdy,!`),n), sublist); 2886 n:=n-1; 2887 xlist:=sqrest xlist 2888 >>; 2889 return sublist 2890end$ 2891 2892%------------- 2893 2894symbolic operator totdeg$ 2895symbolic procedure totdeg(p,f)$ 2896% Ordnung (total) der hoechsten Ableitung von f im Ausdruck p 2897eval(cons('plus,ldiff1(car ldifftot(reval p,reval f),fctargs reval f)))$ 2898 2899%------------- 2900 2901symbolic procedure combidif(s)$ 2902% extracts the list of derivatives from s: % u`1`1`2 --> (u, 1, 1, 2) 2903begin scalar temp,ans,no,n1; 2904 s:=reval s; % to guarantee s is in true prefix form 2905 temp:=reverse explode s; 2906 2907 while not null temp do 2908 <<n1:=<<no:=nil; 2909 while (not null temp) and (not eqcar(temp,'!`)) do 2910 <<no:=car temp . no;temp:=cdr temp>>; 2911 compress no 2912 >>; 2913 if (not fixp n1) then n1:=intern n1; 2914 ans:=n1 . ans; 2915 if eqcar(temp,'!`) then <<temp:=cdr temp; temp:=cdr temp>>; 2916 >>; 2917 return ans 2918end$ 2919 2920%------------- 2921 2922symbolic operator dif$ 2923symbolic procedure dif(s,n)$ 2924% e.g.: dif(fnc!`1!`3!`3!`4, 3) --> fnc!`1!`3!`3!`3!`4 2925begin scalar temp,ans,no,n1,n2,done; 2926 s:=reval s; % to guarantee s is in true prefix form 2927 temp:=reverse explode s; 2928 n2:=reval n; 2929 n2:=explode n2; 2930 2931 while (not null temp) and (not done) do 2932 <<n1:=<<no:=nil; 2933 while (not null temp) and (not eqcar(temp,'!`)) do 2934 <<no:=car temp . no;temp:=cdr temp>>; 2935 compress no 2936 >>; 2937 if (not fixp n1) or ((fixp n1) and (n1 leq n)) then 2938 <<ans:=nconc(n2,ans); ans:='!` . ans; ans:='!! . ans; done:=t>>; 2939 ans:=nconc(no,ans); 2940 if eqcar(temp,'!`) then <<ans:='!` . ans; ans:='!! . ans; 2941 temp:=cdr temp; temp:=cdr temp>>; 2942 >>; 2943 return intern compress nconc(reverse temp,ans); 2944end$ 2945 2946%------------- 2947 2948%symbolic operator totdif$ 2949%symbolic procedure totdif(s,x,n,dylist)$ 2950%% total derivative of s(x,dylist) w.r.t. x which is the n'th variable 2951%begin 2952% scalar tdf,el1,el2; 2953% tdf:=simpdf {s,x}; 2954% <<dylist:=cdr dylist; 2955% while dylist do 2956% <<el1:=cdar dylist;dylist:=cdr dylist; 2957% while el1 do 2958% <<el2:=car el1;el1:=cdr el1; 2959% tdf:=addsq(tdf ,multsq( simp!* dif(el2,n), simpdf {s,el2})) 2960% >> 2961% >> 2962% >>; 2963% return prepsq tdf 2964%end$ 2965 2966put('totdif,'psopfn,'tot!*dif)$ 2967 2968symbolic procedure tot!*dif(inp)$ 2969% total derivative of s(x,dylist) w.r.t. x which is the n'th variable 2970begin 2971 scalar tdf,el1,el2,s,x,n,dylist; 2972 s := aeval car inp$ s:=if pairp s and (car s='!*sq) then cadr s 2973 else simp s; 2974 x := reval cadr inp$ if pairp x and (car x='!*sq) then x:=reval x; 2975 n := reval caddr inp$ 2976 dylist:=cdr reval cadddr inp$ 2977 tdf:=diffsq(s,x); 2978 while dylist do 2979 <<el1:=cdar dylist;dylist:=cdr dylist; 2980 while el1 do 2981 <<el2:=car el1;el1:=cdr el1; 2982 tdf:= addsq(tdf, multsq( simp!* dif(el2,n), diffsq(s,el2))) 2983 >> 2984 >>; 2985 return {'!*sq, tdf, t} 2986end$ 2987 2988endmodule$ 2989 2990%******************************************************************** 2991module divintegration$ 2992%******************************************************************** 2993% Routines to write a given expression as divergence 2994% Author: Thomas Wolf 2995% 1998 2996 2997%------------- 2998 2999%symbolic operator combi$ 3000symbolic procedure combi(ilist)$ 3001% ilist is a list of indexes (of variables of a partial derivative) 3002% and returns length!/k1!/k2!../ki! where kj! is the multiplicity of j. 3003begin 3004 integer n0,n1,n2,n3; 3005 n1:=1; 3006% ilist:=cdr ilist; 3007 while ilist do 3008 <<n0:=n0+1;n1:=n1*n0; 3009 if car ilist = n2 then <<n3:=n3+1; n1:=n1/n3>> 3010 else <<n2:=car ilist; n3:=1>>; 3011 ilist:=cdr ilist>>; 3012 return n1 3013end$ 3014 3015%------------- 3016 3017symbolic procedure sortli(l)$ 3018% sort a list of numbers 3019begin scalar l1,l2,l3,m,n$ 3020 return 3021 if null l then nil 3022 else << 3023 n:=car l$ 3024 l2:=list car l$ 3025 l:=cdr l$ 3026 while l do << 3027 m:=car l$ 3028 if m>n then l1:=cons(car l,l1) 3029 else if m<n then l3:=cons(car l,l3) 3030 else l2:=cons(car l,l2)$ 3031 l:=cdr l 3032 >>$ 3033 append(sortli(l1),append(l2,sortli(l3))) 3034 >> 3035end$ 3036 3037%------------- 3038 3039symbolic procedure derili(il)$ 3040% make a derivative index list from a list of numbers 3041if null il then nil else 3042begin scalar h1,h2,h3$ 3043 h1:=sortli(il); 3044 while h1 do << 3045 h2:=reval algebraic mkid(!`,lisp car h1); 3046 h3:=if h3 then mkid(h2,h3) 3047 else h2; 3048 h1:=cdr h1 3049 >>; 3050 return h3 3051end$ 3052 3053%------------- 3054 3055symbolic procedure newil(il,mo,nx)$ 3056if (null il) or (length il<mo) then cons(1,il) else 3057if car il<nx then cons(add1 car il,cdr il) else 3058<<while il and (car il = nx) do il:=cdr il; 3059 if null il then nil 3060 else cons(add1 car il,cdr il)>>$ 3061 3062%------------- 3063 3064symbolic operator intcurrent1$ 3065% used in conlaw2,4 3066symbolic procedure intcurrent1(divp,ulist,xlist,dulist, 3067 nx,eqord,densord)$ 3068% computes a list in reverse order from which the conserved 3069% current is computed through integration 3070begin scalar ii,il,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11,contrace,u, 3071 nega,pii,mo,pl,nu$ 3072 % contrace:=t; 3073 3074 xlist:=cdr xlist; 3075 ulist:=cdr ulist; 3076 nu:=length ulist; 3077 mo:=if eqord>densord then eqord-1 3078 else densord-1$ 3079 3080 pl:=for ii:=1:nx collect << % the components W^i 3081 il:=nil; 3082 pii:=nil; 3083 repeat << 3084 h11:=cons(ii,il); 3085 h1:=derili(h11); 3086 h11:=combi(sortli(h11)); 3087 if contrace then 3088 <<write"==== ii=",ii," il=",il," h1=",h1," h11=",h11;terpri()>>; 3089 h2:=il; 3090 h3:=nil; 3091 if null h2 then h4:=list(nil . nil) 3092 else << 3093 h4:=list(car h2 . nil); 3094 while h2 do << 3095 h3:=cons(car h2,h3);h2:=cdr h2; 3096 h4:=cons((if h2 then car h2 3097 else nil ) . derili(h3),h4)$ 3098 >>; 3099 >>; 3100 if contrace then <<write"h3=",h3," h4=",h4;terpri()>>; 3101 for e1:=1:nu do << % for each function u 3102 u:=nth(ulist,e1); 3103 h5:=mkid(u,h1); 3104 if contrace then <<write"h5=",h5;terpri()>>; 3105% h6:=nil; 3106 if contrace then <<write"h6-1=", list('df,divp,h5); 3107 terpri()>>; 3108% h6:=cons(reval list('quotient,list('df,divp,h5),h11), h6); 3109% if nde=1 then h6:=car h6 3110% else h6:=cons('plus,h6); 3111 h6:=reval list('quotient,list('df,divp,h5),h11); 3112 if contrace then <<write"h6=",h6;terpri()>>; 3113 nega:=nil; 3114 h7:=h4; 3115 % h7 is in such an order that the first term is always positive 3116 3117 while h7 and not zerop h6 do << 3118 h8:=car h7; h7:=cdr h7; 3119 h9:=car h8; h8:=cdr h8; 3120 if contrace then <<write if nega then "--" else "++", 3121 " h8=",h8," h9=",h9;terpri()>>; 3122 if contrace then <<write"h9-2=",h9;terpri(); 3123 write"h6=",h6$terpri(); 3124 write"xlist=",xlist$terpri(); 3125 write"dulist=",dulist$terpri()>>; 3126 if h9 then h6:=reval tot!*dif({aeval h6,nth(xlist,h9),h9,dulist}); 3127 3128 if contrace then <<write"h6-3=",h6;terpri()>>; 3129 3130 h10:=if h8 then mkid(u,h8) 3131 else u; 3132 if contrace then <<write"h10=",h10;terpri()>>; 3133 h10:=list('times,h6,h10); 3134 if nega then h10:=list('minus,h10); 3135 if contrace then algebraic write">>>> h10=",h10; 3136 pii:=cons(h10,pii)$ 3137 nega:=not nega; 3138 >> 3139 >>; % for each function u 3140 il:=newil(il,mo,nx); 3141 >> until null il; 3142 pii:=reval if null pii then 0 else 3143 if length pii=1 then car pii 3144 else cons('plus,pii); 3145 if contrace then algebraic write"pii-end=",pii; 3146 3147 pii 3148 >>; % for all ii 3149 return cons('list,pl) 3150end$ % of intcurrent1 3151 3152%------------- 3153 3154symbolic operator intcurrent2$ 3155% used in conlaw2,4, crdec 3156symbolic procedure intcurrent2(divp,ulist,xlist)$ 3157% computes the conserved current P_i from the divergence through 3158% partial integration 3159% potential improvement: one could substitute low order derivatives 3160% by high order derivatives using remaining conditions and try again 3161 3162begin scalar h2,h3,h4,h5,h6,h7,h8,e2,e3,lin_problem_bak; 3163 lin_problem_bak:=lin_problem$ lin_problem:=t; % for intpde 3164 3165 % repeated partial integration to compute P_i 3166 ulist :=cdr reval ulist; 3167 xlist :=cdr reval xlist; 3168 3169 h4:=list xlist; 3170 % dequ is here a list containing only the conserved density 3171 % and flux to drop commen factors 3172 repeat << 3173 e3:=divp; 3174 h3:=car h4; % h3 is the list of variables is a spec. order 3175 h4:=cdr h4; 3176 h5:=for e2:=1:length h3 collect 0; 3177 % h5 is old list of the conserved quantities 3178 h8:=0; % h8 counts integrations wrt. all variables 3179 repeat << % integrate several times wrt. all variables 3180 h8:=h8+1; 3181 h2:=h3; % h2 is a copy of h3 3182 h6:=nil; % h6 is new list of the conserved quantities 3183 h7:=nil; % h7 is true if any integration was possible 3184 while h2 neq nil do << % integrating wrt. each variable 3185 e2:=intpde(e3,ulist,h3,car h2,t); 3186 if null e2 then e2:=list(nil,e3) 3187 else e3:=cadr e2; 3188 if (car e2) and (not zerop car e2) then h7:=t; 3189 h6:=cons(list('plus,car e2,car h5),h6); 3190 h5:=cdr h5; 3191 h2:=cdr h2 3192 >>; 3193 h5:=reverse h6; 3194 >> until (h7=nil) % no integration wrt. no variable was possible 3195 or (e3=0) % complete integration 3196 or (h8=10); % max. 10 integrations wrt. all variables 3197 >> until (e3=0) or (h4=nil); 3198 lin_problem:=lin_problem_bak; 3199 3200 return {'list,reval cons('list, h5),e3} 3201 % end of the computation of the conserved current P 3202 % result is a {{P_i},remainder} 3203 % was successful if 0 = remainder (=e3) 3204end$ % of intcurrent2 3205 3206 3207%------------- 3208 3209symbolic operator intcurrent3$ 3210% crident 3211symbolic procedure intcurrent3(divp,ulist,xlist)$ 3212% computes the conserved current P_i from the divergence through 3213% partial integration with restriction of maximal 2 terms 3214 3215begin scalar xl,h1,h2,h3,h4,h5,resu1,resu2,resu3,succ,lin_problem_bak; 3216 lin_problem_bak:=lin_problem$ lin_problem:=t; % for intpde 3217 3218 % repeated partial integration to compute P_i 3219 ulist :=cdr reval ulist; 3220 xlist :=cdr reval xlist; 3221 xl:=xlist; 3222 resu1:=nil; 3223 succ:=nil; 3224 % try all possible different pairs of variables 3225 while (cdr xl) and not succ do << 3226 h1:=intpde(divp,ulist,xlist,car xl,t); 3227 if h1 and not zerop car h1 then << 3228 resu2:=cons(car h1,resu1); 3229 h2:=cdr xl; 3230 repeat << 3231 h3:=intpde(cadr h1,ulist,xlist,car h2,nil); 3232 if h3 and zerop cadr h3 then << 3233 h4:=cons(car h3,resu2); 3234 for each h5 in cdr h2 do h4:=cons(0,h4); 3235 succ:=t; 3236 resu3:= {'list,cons('list,reverse h4),0} 3237 >>; 3238 h2:=cdr h2; 3239 resu2:=cons(0,resu2) 3240 >> until succ or null h2; 3241 >>; 3242 resu1:=cons(0,resu1); 3243 xl:=cdr xl 3244 >>$ 3245 3246 lin_problem:=lin_problem_bak; 3247 return if succ then resu3 3248 else {'list,cons('list,cons(0,resu1)),divp} 3249end$ % of intcurrent3 3250 3251endmodule$ 3252 3253%******************************************************************** 3254module quasilinpde_integration$ 3255%******************************************************************** 3256% Routines to solve a quasilinear first order PDE 3257% Author: Thomas Wolf 3258% summer 1995 3259 3260%---------------------------- 3261 3262algebraic procedure select_indep_var(clist,xlist)$ 3263begin scalar s,scal,notfound,cq,x,xx,xanz,sanz,xcop,ok; 3264 % Find the simplest non-zero element of clist 3265 notfound:=t; 3266 xcop:=xlist; 3267 while notfound and (xcop neq {}) do << 3268 x :=first xcop ; xcop :=rest xcop ; 3269 cq:=first clist; clist:=rest clist; 3270 if cq neq 0 then << 3271 xanz:=0; 3272 for each xx in xlist do 3273 if %df(cq,xx) 3274 not freeof(cq,xx) then xanz:=xanz+1; 3275 3276 ok:=nil; 3277 if not s then ok:=t else % to guarantee s neq nil 3278 if xanz<sanz then ok:=t else % as few dependencies as poss. 3279 if xanz=sanz then 3280 if length(cq)=1 then ok:=t else 3281 if length(scal)>1 then 3282 if part(scal,0)=plus then % if possible not a sum 3283 if (part(cq,0) neq plus) or 3284 (length(cq)<length(scal)) then ok:=t; 3285 3286 if ok then <<s:=x; scal:=cq; sanz:=xanz>>; 3287 if scal=1 then notfound:=nil 3288 >> 3289 >>; 3290 3291 return {s,scal} 3292end$ % of select_indep_var 3293 3294%---------------------------- 3295 3296algebraic procedure charsyscrack(xlist,cqlist,s,scal,u,ode)$ 3297% formulation and solution of the characteristic ODE-system 3298begin 3299 scalar lcopy,x,cS,flist,soln,printold,timeold,facintold, 3300 adjustold,safeintold,freeintold,odesolveold, 3301 e1,e2,e3,e4,n,nmax,dff,proclistold,flin_old, 3302 lin_problem_old,session_old,level_old; 3303 % formulating the ode/pde - system . 3304 lcopy := cqlist ; 3305 cS := {} ; 3306 for each x in xlist do 3307 << cq := first lcopy ; lcopy := rest lcopy ; 3308 if x neq s then 3309 << depend x,s; 3310 cS := .(scal*df(x,s)-cq,cS) >> 3311 >>; 3312 if s neq u then <<flist:={u}; depend u,s; 3313 cS := .(scal*df(u,s)+ode,cS) >> 3314 else flist:={}; 3315 for each x in xlist do 3316 if s neq x then flist:=.(x,flist); 3317 3318 lisp << 3319 timeold := time_; time_ :=nil; 3320 facintold:=facint_; facint_:=1000; 3321 adjustold:=adjust_fnc; adjust_fnc:=t; 3322 safeintold:=safeint_; safeint_:=nil; 3323 freeintold:=freeint_; freeint_:=nil; 3324 odesolveold:=odesolve_; odesolve_:=50; 3325 proclistold:=proc_list_; 3326 proc_list_:=delete('stop_batch,delete('alg_solve_single,proc_list_))$ 3327 flin_old:=flin_$ 3328 lin_problem_old:=lin_problem$ 3329 session_old:=session_$ 3330 level_:=level_old; 3331 >>$ 3332 % solving the odes using crack. 3333 if lisp(print_ neq nil) then 3334 lisp << 3335 write "The equivalent characteristic system:";terpri(); 3336 deprint cdr algebraic cS$ 3337 %terpri()$ 3338 write "for the functions: "; 3339 fctprint( cdr reval algebraic flist);write"."; 3340% write"---------------------------------------------------------------------------"$terpri()$ 3341 write"A new CRACK computation will start now to solve a characteristic ODE system."$terpri()$ 3342 >>; 3343 soln:=crack(cS,flist,flist,{}); 3344 lcopy:={}; 3345 for each x in soln do << 3346 e1:=first x; 3347 if (e1 neq {}) and (length e1 + length second x = length flist) then << 3348 3349 % are there enough constants of integration? 3350 e2:=length third x; 3351 for each e3 in flist do 3352 if not freeof(third x,e3) then e2:=e2-1; 3353 if e2=length cS then << % enough constants 3354 3355 % all remaining conditions are algebraic (non-differential) in all the 3356 % functions of flist? 3357 e2:={}; 3358 e3:={}; 3359 for each e4 in third x do if freeof(flist,e4) then e3:=e4 . e3 3360 else e2:=e4 . e2; 3361 if (length cS) = (length e3) then << % sufficiently many integrations done 3362 for each e4 in e2 do 3363 for each e5 in e1 do 3364 if lisp(not freeof(lderiv(reval algebraic e5, 3365 reval algebraic e4, 3366 list(reval algebraic s)),'df)) 3367 then <<e2:={};e1:={}>>; 3368 3369 if e2 neq {} then << 3370 % It may be possible that derivatives of unsolved functions 3371 % occur in the expressions of the evaluated functions: second x 3372 nmax:=0; 3373 for each e4 in e2 do << % for each unsolved function 3374 for each e5 in second x do << % for each solved expression 3375 lisp << 3376 n:=lderiv(reval algebraic rhs e5,reval algebraic e4, 3377 list(reval algebraic s)); 3378 n:=if (not pairp car n) or (caar n neq 'df) then 0 else 3379 % would be wrong, for example, for n=sin(df(..,..)) 3380 if (length car n) = 3 then 1 3381 else cadddr car n 3382 >>; 3383 n:=lisp n; 3384 if n>nmax then nmax:=n; 3385 >>; 3386 if nmax>0 then << % adding further conditions 3387 e5:=e1; 3388 while freeof(first e5,e4) do e5:=rest e5; 3389 e5:=first e5; 3390 dff:=e4; 3391 for n:=1:nmax do << 3392 e5 :=df(e5,s); e1:= e5 . e1; 3393 dff:=df(dff,s); e3:=dff . e3 3394 >> 3395 >> 3396 >>; 3397 lcopy:=cons({append(second x,e1),e3},lcopy); 3398 >> 3399 >> 3400 >> 3401 >> else 3402 if (first x = {}) and (length cS = length third x) then 3403 lcopy:=cons({second x,third x},lcopy) 3404 >>; 3405 3406 lisp << 3407 time_:=timeold; 3408 facint_:=facintold; 3409 adjust_fnc:=adjustold; 3410 safeint_:=safeintold; 3411 freeint_:=freeintold; 3412 odesolve_:=odesolveold; 3413 proc_list_:=proclistold; 3414 flin_:=flin_old; 3415 lin_problem:=lin_problem_old; 3416 session_:=session_old; 3417 level_:=level_old; 3418 >>; 3419 return 3420 if lcopy={} then <<for each x in flist do 3421 if not my_freeof(x,s) then nodepend x,s; {}>> 3422 else s . lcopy 3423 % { {{x=..,y=..,u=..,..,0=..},{df(z,s),df(z,s,2),..,c1,c2,c3,..}},..} 3424 % df(z,s,..) only if df(z,s,..) turns up in x, y, u. .. . 3425end$ % of charsyscrack 3426 3427%---------------------------- 3428 3429%procedure charsyspsfi(xlist,cqlist,u,ode,denf); 3430% % This calls psfi(). Not sure where this is defined (16.12.2007) 3431%begin 3432% scalar h,s; 3433% h:=cqlist; 3434% cqlist:={}; 3435% while h neq {} do <<cqlist:=cons(first h*denf,cqlist);h:=rest h>>; 3436% cqlist:=cons(-ode*denf,cqlist); 3437% xlist:=cons(u,reverse xlist); 3438% s:=lisp gensym(); 3439% for each h in xlist do depend h,s; 3440% h:=psfi(cqlist,xlist,s); 3441% for each h in xlist do if not my_freeof(h,s) then nodepend h,s; 3442% return h 3443%end$ % of charsyspsfi 3444 3445%---------------------------- 3446 3447algebraic procedure storedepend(xlist)$ 3448% stores the dependencies of the elements of xlist in the returned 3449% list and clears the dependencies 3450begin scalar q,e1,e2$ 3451 return 3452 for each e1 in xlist collect 3453 << q:=fargs e1; 3454 for each e2 in q do nodepend e1,e2; 3455 cons(e1,q)>> 3456end$ % of storedepend 3457 3458%---------------------------- 3459 3460algebraic procedure restoredepend(prev_depend)$ 3461% assigns the dependencies stored in the argument 3462begin scalar q,s,x; 3463 for each s in prev_depend do 3464 <<q:=first s; s:=rest s; 3465 for each x in s do depend q,x>> 3466end$ % of restoredepend 3467 3468%---------------------------- 3469 3470symbolic procedure dropable(h,fl,allowed_to_drop)$ 3471if (not smemberl(fl,h)) or 3472 member(h,allowed_to_drop) then t else nil$ 3473 3474%---------------------------- 3475 3476symbolic procedure drop_factors(q,fl,allowed_to_drop)$ 3477if not pairp q or (car q neq 'times) then 3478if dropable(q,fl,allowed_to_drop) then 1 3479 else q else 3480reval cons('times,for each h in cdr q collect 3481 if dropable(h,fl,allowed_to_drop) then 1 else h)$ 3482 3483%---------------------------- 3484 3485symbolic procedure simplifiziere(q,fl,allowed_to_drop)$ 3486begin 3487 scalar n,nu,de; 3488 return 3489 if not pairp q then q else 3490 if member(car q,ONE_ARGUMENT_FUNCTIONS_) or 3491 (member(car q,{'expt}) and 3492 dropable(caddr q,fl,allowed_to_drop)) then simplifiziere(cadr q,fl,allowed_to_drop) 3493 else 3494 if car q = 'quotient then << 3495 nu:=drop_factors(cadr q,fl,allowed_to_drop); 3496 de:=drop_factors(caddr q,fl,allowed_to_drop); 3497 if nu=1 then simplifiziere(de,fl,allowed_to_drop) else 3498 if de=1 then simplifiziere(nu,fl,allowed_to_drop) else {'quotient,nu,de} 3499 >> else 3500 if car q = 'times then << 3501 de:=drop_factors(q,fl,allowed_to_drop); 3502 if (not pairp de) or (length de < length q) then 3503 simplifiziere(de,fl,allowed_to_drop) else de 3504 >> else << 3505 q:=signchange(q); 3506 n:=ncontent(q); 3507 if n=1 then q % one could also subtract constant multiples of 3508 % elements of allowed_to_drop to simplify q 3509 else simplifiziere(reval list('quotient, q, reval n), 3510 fl,allowed_to_drop) 3511 >> 3512end$ 3513 3514%---------------------------- 3515 3516algebraic procedure quasilinpde(f,u,xlist)$ 3517% f ... PDE, u ... unknown function, xlist ... variable list 3518begin scalar i, q, qlist, cq, cqlist, ode, soln, tran, 3519 xcop, s, s1, x, xx, h1, h2, scal, qlin, prev_depend, 3520 xlist_cp1,xlist_cp2; 3521 symbolic put('ff,'simpfn,'simpiden)$ 3522 if lisp print_more then 3523 write"The quasilinear PDE: 0 = ",f,"."; 3524 % changing the given pde into a quasi-linear ode . 3525 i := 0 ; 3526 ode := f; 3527 qlist := {}; cqlist:={}; 3528 for each x in xlist do 3529 <<i := i+1 ; 3530 q := mkid(p_,i) ; 3531 qlist := .(q,qlist) ; 3532 f := sub(df(u,x) = q , f) ; 3533 cq:=df(f,q); 3534 cqlist:=.(df(f,q),cqlist) ; 3535 ode := ode - df(u,x)*df(f,q) ; 3536 %if not my_freeof(u,x) then nodepend u, x ; 3537 >> ; 3538 lisp(depl!*:=delete(assoc(u,depl!*),depl!*))$ 3539 lisp(depl!*:=delete(assoc(mkid(u,'_),depl!*),depl!*))$ 3540 3541 qlist := reverse qlist ; 3542 cqlist := reverse cqlist ; 3543 3544 % checking for linearity. 3545 qlin:=t; 3546 for each cq in cqlist do 3547 for each q in qlist do 3548 if df(cq,q) neq 0 then qlin:=nil; 3549 if not qlin then return {}; 3550 3551% soln:=charsyspsfi(xlist,cqlist,u,ode,den f); 3552 3553 % Determination of the independent variable for the ODEs 3554 pcopy:=cons(-ode,cqlist);xcop:=cons(u,xlist); 3555 scal:=select_indep_var(pcopy,xcop)$ 3556 s1:=first scal; 3557 3558 prev_depend:=storedepend(xlist)$ 3559 soln:=charsyscrack(xlist,cqlist,s1,second scal,u,ode); 3560 if soln={} then << % try all other coefficients as factors 3561 repeat << 3562 repeat <<s :=first xcop ;xcop :=rest xcop; 3563 scal:=first pcopy;pcopy:=rest pcopy>> 3564 until (pcopy={}) or ((scal neq 0) and (s neq s1)); 3565 if (s neq s1) and (scal neq 0) then << 3566 if lisp print_ then lisp 3567 <<terpri()$ 3568 write"New attempt with a different independent variable:"; 3569 terpri()>>$ 3570 soln:=charsyscrack(xlist,cqlist,s,scal,u,ode) 3571 >> 3572 >> 3573 until (soln neq {}) or (xcop={}) 3574 >>; 3575 % solving for the constants(eg..c1,c2 etc.) and put them in a 3576 % linear form. 3577 % The first element of soln is the ODE-variable 3578 tran:={}; 3579 if soln neq {} then << 3580 s1:=first soln; 3581 for each s in rest soln do 3582 <<cq:=first solve(first s,second s); 3583 x:={}; 3584 for each xx in cq do 3585 if lisp atom algebraic lhs xx then 3586 x:=.(lisp <<q:=reval algebraic rhs xx; 3587 simplifiziere(q,cons(reval u,cdr xlist),nil)>>,x); 3588 lisp << 3589 x:=cdr x; 3590 repeat << 3591 xx:=x$ 3592 x:=for each h1 in x collect 3593 simplifiziere(h1,cons(reval u,cdr xlist),delete(h1,x)) 3594 >> until x=xx$ 3595 x:=cons('list,x); 3596 >>$ 3597 xx:=tran; 3598 while (xx neq {}) and 3599 <<h1:=x;h2:=first xx; 3600 while (h1 neq {}) and 3601 (h2 neq {}) and 3602 (first h1 = first h2) do <<h1:=rest h1; h2:=rest h2>>; 3603 if (h1={}) and (h2={}) then nil 3604 else t 3605 >> do xx:=rest xx; 3606 if xx={} then tran:=.(x,tran); 3607 >>; 3608 for each s in xlist do 3609 if (s neq s1) and (not my_freeof(s,s1)) then nodepend s,s1 3610 >>; 3611 3612 for each x in xlist do depend u,x; 3613 3614 if lisp print_more then 3615 if tran neq {} then << 3616 write"The general solution of the PDE is given through"; 3617 for each x in tran do write"0 = ", 3618 lisp( cons('ff,cdr reval algebraic x)); 3619 if length tran>1 then write"with arbitrary function(s) ff(..)." 3620 else write"with arbitrary function ff(..)." 3621 >>; 3622 % restoring the dependencies of the variables of the PDE 3623 restoredepend(prev_depend)$ 3624 3625 return tran 3626end$ % of quasilinpde 3627 3628endmodule$ 3629 3630end$ 3631 3632integration 3633 integrate_one_pde 3634 integrate 3635 integratede 3636 integrableode 3637 integrableode1 3638 odecheck 3639 odetest 3640 integrateode 3641 odeconvert 3642 odesolve 3643 integratepde 3644 addintco 3645 ld_deriv_search 3646 potintegrable 3647 multipleint 3648 intpde 3649 intpde_ 3650 3651tr integration 3652tr integrate_one_pde 3653tr integrate 3654tr integratede 3655tr integrableode 3656tr integrableode1 3657tr odecheck 3658tr odetest 3659tr integrateode 3660tr odeconvert 3661tr err_catch_odesolve 3662tr odesolve 3663tr integratepde 3664tr addintco 3665tr ld_deriv_search 3666tr potintegrable 3667tr multipleint 3668tr intpde 3669tr intpde_ 3670 367130.1.2015: 3672 3673 - It can not be that case conditions are added which involve new functions 3674 which are not passed on and which are not available after the end of this 3675 integration module. 3676 --> check all places 3677 3678 - It must not happen, that a sub-case is started which does not change the 3679 equation from which the case distinction resulted and then integration is 3680 tried again and the same case distinction with new and different constants 3681 of integration is started again and again. What is needed is that before 3682 case distinctions are written into to_do_list, to check that these 3683 conditions are not already satisfied. 3684 3685 - A fix could be to return the integrated equation multiplied by all cases 3686 Which changes are needed for that? 3687