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