1module cde_diffcon; % CDE package, computation of the
2                    % differential consequences of the initial system of PDEs.
3
4% Redistribution and use in source and binary forms, with or without
5% modification, are permitted provided that the following conditions are met:
6%
7%    * Redistributions of source code must retain the relevant copyright
8%      notice, this list of conditions and the following disclaimer.
9%    * Redistributions in binary form must reproduce the above copyright
10%      notice, this list of conditions and the following disclaimer in the
11%      documentation and/or other materials provided with the distribution.
12%
13% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
14% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
15% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
16% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
17% CONTRIBUTORS
18% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
19% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
20% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
21% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
22% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
23% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
24% POSSIBILITY OF SUCH DAMAGE.
25%
26
27% % *****************************************************************
28% Author and maintainer: Raffaele Vitolo
29% Dipartimento di Matematica, Universita' del Salento (Lecce, Italy)
30% email: raffaele.vitolo@unisalento.it
31% web: http://poincare.unisalento.it/vitolo
32% ===============================================================
33
34%-----------------------------------------------------------------------------%
35% Computation of differential consequences of E=0
36%-----------------------------------------------------------------------------%
37%
38% We compute differential consequences as total derivatives
39% (on the free jetspace!) of de and de_odd,
40% then we solve the system of differential consequences
41% in terms of parametric derivatives only.
42% While solving, other principal derivatives could be needed and computed.
43
44symbolic procedure order_diffcon(par,dcon_temp);
45   % Given an alist of differential consequences it will be ordered
46   % according with decreasing order of principal derivative.
47   % It uses quicksort.
48   % The only purpose of this function is to order initial equations
49   % according with decreasing order of derivatives.
50begin
51  scalar locless,locequal,locgreater,locpivot;
52  if not(dcon_temp) then return nil else return
53  <<
54    locpivot:=car dcon_temp;
55    for each el in dcon_temp do
56      if lessp(order_of_der(par,car el),order_of_der(par,car locpivot))
57        then locless:=el . locless
58      else if eqn(order_of_der(par,car el),order_of_der(par,car locpivot))
59        then locequal:=el . locequal
60      else locgreater:=el . locgreater;
61    append(
62      append(locequal,order_diffcon(par,locgreater)),
63      order_diffcon(par,locless)
64    	)
65  >>;
66end;
67
68symbolic procedure diff_vars(var1,var2);
69   % takes the difference of two variables in multiindex notation;
70   % returns -1 if the dependent part is different
71   % or if the difference of multiindexes contains negative entries;
72   % return the difference of the corresponding multiindexes in other cases.
73  begin scalar m_n,mmn; integer mmni;
74    if neq(car var1,car var2) then return (-1);
75    m_n:=pair(cadr var1,cadr var2);
76    for each el in m_n do
77    <<
78      mmni:=car el - cdr el;
79      if mmni >=0 then mmn:=mmni . mmn
80    >>;
81    if length(mmn)<n_indep_var then return (-1)
82    else return reverse mmn
83  end;
84
85symbolic procedure add_diffcon(par,dcon);
86   % Adds the expression of the differential consequences
87   % corresponding to the list of principal derivatives `dcon'
88   % We assume that diffcon_der!* and diffcon_odd!*
89   % are made of derivatives in decreasing order
90   begin scalar diffcontemp,vartemp,pder,exprtemp; integer ordtemp,dmind,cnt;
91      if par=0 then diffcontemp:=diffcon_der!* else diffcontemp:=diffcon_odd!*;
92      for each el in dcon do
93      <<
94         vartemp:=idtomind(par,el);
95         ordtemp:=order_of_der(par,el);
96         dmind:=-1;
97         cnt:=1;
98         % Next loop searches the starting point of differential consequence
99         % of order less than the order of the variable `el'.
100         while order_of_der(par,car nth(diffcontemp,cnt))>=ordtemp do
101            cnt:=cnt+1;
102         % Next loop searches for a variable with the same dependent part
103         % of the given one and a multiindex which has positive difference
104         % with the given one.
105         while dmind=-1 do
106         <<
107	   pder:=car nth(diffcontemp,cnt);
108           dmind:=diff_vars(vartemp,idtomind(par,pder));
109           cnt:=cnt+1
110         >>;
111         exprtemp:=cadr nth(diffcontemp,cnt-1);
112         for i:=1:n_indep_var do
113            for j:=1:nth(dmind,i) do
114	      exprtemp:=aeval list(nth(tot_der!*,i),exprtemp);
115         if not freeof(exprtemp,'letop) then exprtemp:=aeval 'letop;
116	 % we insert the differential consequence
117         % in such a way to preserve the ordering
118         cnt:=1;
119         while order_of_der(par,el) <
120	    order_of_der(par,car nth(diffcontemp,cnt)) do cnt:=cnt+1;
121         diffcontemp:=cde_insert(list(el,exprtemp),diffcontemp,cnt);
122      >>;
123      if par=0 then diffcon_der!*:=diffcontemp
124      else diffcon_odd!*:=diffcontemp;
125      if par=0 then
126         diffcon_comp_der!*:=for each el in diffcon_der!* collect car el
127      else
128         diffcon_comp_odd!*:=for each el in diffcon_odd!* collect car el;
129   end;
130
131symbolic procedure check_diffcon(dc_der,dc_odd);
132  begin scalar dc_extra_der,dc_extra_odd,dc_new_der,dc_new_odd,
133      ell_ext,exprtemp;
134    % check if computed differential consequences
135    % depend on non-computed (extra) differential consequences.
136    % If such non-computed differential consequences exist, compute them.
137    dc_extra_der:=cde_diffset(all_principal_der!*,diffcon_comp_der!*);
138    dc_extra_odd:=cde_diffset(all_principal_odd!*,diffcon_comp_odd!*);
139    for each el in dc_der do
140    <<
141      exprtemp:=cadr assoc(el,diffcon_der!*);
142      for each ell in dc_extra_der do
143	if not freeof(exprtemp,ell) then dc_new_der:=ell . dc_new_der
144    >>;
145    for each el in dc_odd do
146    <<
147      exprtemp:=cadr assoc(el,diffcon_odd!*);
148      for each ell in dc_extra_der do
149	if not freeof(exprtemp,ell) then dc_new_der:=ell . dc_new_der;
150      for each ell in dc_extra_odd do
151      <<
152	ell_ext:=oddext(ell);
153	if not freeof(exprtemp,ell_ext) then
154	  dc_new_odd:=ell . dc_new_odd
155      >>
156    >>;
157    return list(cde_mkset(dc_new_der),cde_mkset(dc_new_odd))
158  end;
159
160symbolic procedure generate_diffcon();
161   % Starting from the primary differential consequences, the code:
162   % 1 adds all expressions of differential consequences;
163   % 2 makes a scan of them in order to find
164   %   secondary differential consequences;
165   % 3 if there are secondary diff. cons, it adds all differential consequences
166   %   and goes to item 2.
167   begin scalar dcon_der,dcon_odd,tempder; integer dependency,n_dcon;
168      % Initialization of the differential consequences.
169      % These are lists of even (odd) principal derivatives
170      % paired with their expression.
171      % The lists are ordered according with decreasing order of derivatives.
172      diffcon_der!*:=
173	order_diffcon(0,pair(principal_der!*,
174	  for each el in de!* collect list(el)));
175      diffcon_odd!*:=
176	order_diffcon(1,pair(principal_odd!*,
177	  for each el in de_odd!* collect list(replace_oddext(el))));
178      diffcon_comp_der!*:=for each el in diffcon_der!* collect car el;
179      diffcon_comp_odd!*:=for each el in diffcon_odd!* collect car el;
180      dependency:=1;
181      dcon_der:=primary_diffcon_der!*;
182      dcon_odd:=primary_diffcon_odd!*;
183      n_dcon:=1;
184      while dependency=1 do
185      <<
186         % compute current differential consequences;
187         % it could be that one of the current dcon files is empty.
188         if dcon_der then add_diffcon(0,dcon_der);
189         if dcon_odd then add_diffcon(1,dcon_odd);
190         tempder:=check_diffcon(dcon_der,dcon_odd);
191         dcon_der:=car tempder;
192         dcon_odd:=cadr tempder;
193         if and(not(dcon_der),not(dcon_odd)) then dependency:=0
194	 else
195	 <<
196	   n_dcon:=n_dcon+1;
197	   prin2 "Presence of ";prin1 n_dcon;
198	   prin2t "-ary differential consequences";
199	 >>
200      >>;
201   end;
202
203symbolic procedure replace_diffcon_totder(even_diffcon_tot,odd_diffcon_tot,
204    prin_param_der,prin_param_ext);
205    begin
206      for each el in even_diffcon_tot do
207      	set_svf(car el,0,caddr el,
208          cdr assoc(car reverse el,prin_param_der)
209      	    );
210      for each el in odd_diffcon_tot do
211	set_svf(car el,1,caddr el,
212	  cdr assoc(oddext(car reverse el),prin_param_ext)
213      	    );
214    end;
215
216symbolic procedure cde_differential_consequences();
217  begin scalar templdcon,diffcon_comp_der_paramexpr,diffcon_comp_ext_expr,
218      diffcon_comp_ext_evenparamexpr,diffcon_comp_ext_paramexpr;
219    % We start by generating differential consequences.
220    prin2t "   a - Calculating differential consequences ...";
221    generate_diffcon();
222    prin2t "   b - Solving the system of differential consequences ...";
223    % We solve the system of differential consequences in order
224    % to obtain principal derivatives only on the lhs.
225    % We do it in steps in order to limit memory occupation.
226    % 0 - First of all we prepare new lists of variables
227    %     and corresponding expressions:
228    diffcon_comp_ext!*:=for each el in diffcon_comp_odd!* collect oddext(el);
229    diffcon_comp_ext_expr:=for each el in diffcon_odd!* collect cadr el;
230    %
231    % 1 - replaces all computed even principal derivatives in the rhs
232    %     in order to obtain a 'clean' system:
233    %     even principal derivatives=even parametric expressions;
234    repprincparam_der:=for each el in diffcon_der!* collect
235      list('replaceby,car el,cadr el);
236    diffcon_comp_der_paramexpr:=cdr evalwhereexp(
237      list('list . repprincparam_der,'list . diffcon_comp_der!*)
238    	);
239    diffcon_param_der!*:=pair(diffcon_comp_der!*,diffcon_comp_der_paramexpr);
240    % renews repprincparam_der in order to have a 'clean' replacement list:
241    repprincparam_der:=for each el in diffcon_param_der!* collect
242      list('replaceby,car el,cdr el);
243    %
244    % 2 - replaces all even principal derivatives into the prolongation of
245    %     odd equations and prepares the substitution list;
246%%     prin2t "repprincparam_der:=";
247%%     prin1 repprincparam_der;terpri();
248%%     prin2t "diffcon_comp_ext_expr:";
249%%     prin1 diffcon_comp_ext_expr;terpri();
250    diffcon_comp_ext_evenparamexpr:=
251%      for each el in diffcon_comp_ext_expr collect
252%      cdr evalwhereexp(list('list . repprincparam_der,'list . list(el)));
253      cdr evalwhereexp(
254	list('list . repprincparam_der,'list . diffcon_comp_ext_expr)
255      );
256    %
257    % 3 - replaces all computed odd principal derivatives on the rhs
258    %     in order to obtain a 'clean' system:
259    %     odd principal derivatives=odd parametric expressions;
260    templdcon:=pair(diffcon_comp_ext!*,diffcon_comp_ext_evenparamexpr);
261    repprincparam_ext:=for each el in templdcon collect
262      list('replaceby,car el,cdr el);
263    diffcon_comp_ext_paramexpr:=
264      cdr evalwhereexp(
265	list('list . repprincparam_ext,'list . diffcon_comp_ext!*)
266	  );
267    diffcon_param_ext!*:=
268      pair(diffcon_comp_ext!*,diffcon_comp_ext_paramexpr);
269    % Renew repprincparam_ext in order to have a `clean' replacement list
270    repprincparam_ext:=for each el in diffcon_param_ext!* collect
271      list('replaceby,car el,cdr el);
272    % Define an odd version of repprincparam_ext
273    repprincparam_odd:=for each el in diffcon_param_ext!* collect
274      list('replaceby,extodd(car el),replace_extodd(cdr el));
275    % Redefine repprincparam_der and repprincparam_ext to be algebraic lists
276    repprincparam_der:='list . repprincparam_der;
277    repprincparam_odd:='list . repprincparam_odd;
278    % Finally, we replace the differential consequences
279    % in the expression of total derivatives
280    prin2t "   c - Replacing differential cons. in total derivatives ...";
281    replace_diffcon_totder(primary_diffcon_der_tot!*,primary_diffcon_odd_tot!*,
282      diffcon_param_der!*,diffcon_param_ext!*);
283  end;
284
285algebraic procedure restrict_to_equation(fct);
286  begin
287    return (fct where repprincparam_der)
288  end;
289
290endmodule;
291
292end;
293