1module integrator; 2 3 4% Redistribution and use in source and binary forms, with or without 5% modification, are permitted provided that the following conditions 6% are met: 7% 8% * Redistributions of source code must retain the relevant 9% copyright notice, this list of conditions and the following 10% disclaimer. 11% * Redistributions in binary form must reproduce the above 12% copyright notice, this list of conditions and the following 13% disclaimer in the documentation and/or other materials provided 14% with the distribution. 15% 16% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 17% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 18% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 19% A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 20% OWNERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 21% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 22% LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 23% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 24% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 25% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 26% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 27% 28% % ***************************************************************** 29% 30% Authors: P. Gragert, P.H.M. Kersten, G.H.M. Roelofs, G.F. Post 31% University of Twente (Enschede, The Netherlands) 32% 33% Version and Date: Version 1.0, 1992. 34% 35% Maintainer: Raffaele Vitolo 36% Dipartimento di Matematica, Universita' del Salento (Lecce, Italy) 37% email: raffaele.vitolo@unisalento.it 38% web: http://poincare.unisalento.it/vitolo 39% =============================================================== 40 41symbolic$ 42 43% write"Integrator package for REDUCE, $Revision: 1.0 $"$terpri()$ 44 45put('initialize_equations, 'psopfn, 'initialize_equations1)$ 46 47 48 49global '(cur_eq_set!*)$ 50cur_eq_set!*:= 'equ$ 51 52 53fluid '(!*coefficient_check)$ 54!*coefficient_check:=t$ 55flag('(coefficient_check), 'switch)$ 56 57 58fluid '(!*polynomial_check)$ 59!*polynomial_check:=nil$ 60flag('(polynomial_check), 'switch)$ 61 62 63fluid '(!*allow_differentiation)$ 64!*allow_differentiation:=nil$ 65flag('(allow_differentiation), 'switch)$ 66 67 68 69fluid '(listpri_depth!*)$ 70listpri_depth!*:=40$ 71 72 73algebraic$ 74 75 76lisp procedure initialize_equations1 specification_list; 77 begin scalar operator_name,total_used,variable_list, 78 specification,even_used,odd_used, 79 constant_operator,bracketname,function_name,function_list; 80 if length specification_list<5 then 81 rederr("INITIALIZE_EQUATIONS: wrong number of parameters"); 82 if not idp(operator_name:=car specification_list)then 83 rederr("INITIALIZE_EQUATIONS: equations operator must be identifier"); 84 if not fixp(total_used:= 85 reval car(specification_list:=cdr specification_list)) 86 or total_used<0 then 87 rederr("INITIALIZE_EQUATIONS: total number of equations must be positive"); 88 put(operator_name, 'total_used,total_used); 89 variable_list:=reval car( 90 specification_list:=cdr specification_list); 91 if atom variable_list or car variable_list neq 'list then 92 rederr("INITIALIZE_EQUATIONS: variable list must be algebraic list"); 93 put(operator_name, 'variable_list,cdr variable_list); 94 95 specification_list:=cdr specification_list; 96 specification:=car specification_list; 97 98 if atom specification or length specification neq 4 99 or car specification neq 'list 100 or not idp(constant_operator:=cadr specification)or 101 not fixp(even_used:=reval caddr specification)or 102 not fixp(odd_used:=reval cadddr specification) 103 or even_used<0 or odd_used<0 then 104 105 msgpri("INITIALIZE_EQUATIONS: invalid declaration of", 106 specification,nil,nil,t); 107 put(operator_name, 'constant_operator,constant_operator); 108 if(bracketname:=get(constant_operator, 'bracketname))then 109 put(operator_name, 'bracketname,bracketname); 110 111 if get(constant_operator, 'bracketname)then 112 define_used(bracketname,list('list,even_used,odd_used)) 113 else 114 begin 115 put(constant_operator, 'even_used,even_used); 116 put(constant_operator, 'odd_used,odd_used); 117 end; 118 119 for each function_specification in cdr specification_list do 120 begin 121 122 if atom function_specification or length function_specification neq 4 123 or car function_specification neq 'list 124 or not idp(function_name:=cadr function_specification)or 125 not fixp(even_used:=reval caddr function_specification)or 126 not fixp(odd_used:=reval cadddr function_specification) 127 or even_used<0 or odd_used<0 then 128 129 msgpri("INITIALIZE_EQUATIONS: invalid declaration of", 130 function_specification,nil,nil,t); 131 132 if get(function_name, 'bracketname)then 133 define_used(bracketname,list('list,even_used,odd_used)) 134 else 135 begin 136 put(function_name, 'even_used,even_used); 137 put(function_name, 'odd_used,odd_used); 138 end; 139 function_list:=function_name . function_list; 140 end; 141 put(operator_name, 'function_list,function_list); 142 end$ 143 144 145lisp operator use_equations; 146lisp procedure use_equations operator_name; 147 begin 148 if idp operator_name then 149 cur_eq_set!*:=operator_name 150 else rederr("USE_EQUATIONS: argument must be identifier"); 151 end$ 152 153 154lisp operator integrate_equation; 155lisp procedure integrate_equation n; 156 begin scalar listpri_depth!*,total_used,equation,denominator, 157 solvable_kernel,solvable_kernels,df_list,function_list,present_functions_list, 158 variable_list,absent_variables, 159 linear_functions_list,constants_list,bracketname,df_terms,df_functions, 160 linear_functions,functions_and_constants_list,commutator_functions, 161 present_variables,nr_of_variables,integration_variables; 162 listpri_depth!*:=200; 163 terpri!* t; 164 165 if null(total_used:=get(cur_eq_set!*, 'total_used))or 166 n>total_used then 167 168 msgpri("INTEGRATE_EQUATIONS: properly initialize", 169 cur_eq_set!*,nil,nil,t); 170 if null(equation:=cadr assoc(list(cur_eq_set!*,n), 171 get(cur_eq_set!*, 'kvalue)))then 172 173 msgpri("INTEGRATE_EQUATION:",list(cur_eq_set!*,n), 174 "is non-existent",nil,t); 175 denominator:=denr(equation:=simp!* equation); 176 equation:=numr equation; 177 if null equation then 178 <<write cur_eq_set!*,"(",n,") = 0";terpri!* t; 179 180 setk(list(cur_eq_set!*,n),0);goto solved>> ; 181 182 df_list:=split_form(equation, '(df)); 183 if try_a_homogeneous_integration(n,denominator,df_list)then goto solved; 184 185 186 function_list:=get(cur_eq_set!*, 'function_list); 187 present_functions_list:=get_recursive_kernels(equation,function_list); 188 variable_list:=get(cur_eq_set!*, 'variable_list); 189 absent_variables:=variable_list; 190 for each function in present_functions_list do 191 for each variable in 192 ((if depl_entry then cdr depl_entry) 193 where depl_entry=assoc(function,depl!*))do 194 absent_variables:=delete(variable,absent_variables); 195 if split_equation_polynomially(n,total_used,equation,absent_variables)then 196 goto solved; 197 198 linear_functions_list:=split_form(car df_list, 199 function_list); 200 df_list:=cdr df_list; 201 constants_list:=split_form(car linear_functions_list, 202 list get(cur_eq_set!*, 'constant_operator)); 203 linear_functions_list:=cdr linear_functions_list; 204 if(bracketname:=get(cur_eq_set!*, 'bracketname))then 205 206 if length(df_list)=0 and 207 length(linear_functions_list)=0 then 208 << 209 if atom(solvable_kernel:= 210 relation_analysis(!*ff2a(equation,denominator),bracketname)) 211 then <<write cur_eq_set!*,"(",n,") is a non-solvable Lie relation"; 212 terpri!* t>> 213 else <<write cur_eq_set!*,"(",n,") solved for ";maprin solvable_kernel; 214 terpri!* t; 215 setk(list(cur_eq_set!*,n),0)>> ; 216 goto solved 217 >> ; 218 219 220 df_terms:=for each df_term in df_list join 221 if member(car cadr car df_term,function_list) 222 then list car df_term; 223 for each df_term in df_terms do if not member(cadr 224 df_term,df_functions)then df_functions:=cadr(df_term) . df_functions; 225 functions_and_constants_list:=append(linear_functions_list, 226 cdr constants_list); 227 linear_functions:=for each linear_function in 228 functions_and_constants_list collect car linear_function; 229 if bracketname then commutator_functions:= 230 get_recursive_kernels(car constants_list, 231 get(cur_eq_set!*, 'function_list));; 232 233 present_variables:=variable_list; 234 for each variable in absent_variables do 235 present_variables:=delete(variable,present_variables); 236 nr_of_variables:=length present_variables; 237 for each kernel in linear_functions do if length 238 239 ((if depl_entry then cdr depl_entry) 240 where depl_entry=assoc(kernel,depl!*))=nr_of_variables then 241 solvable_kernels:=kernel . solvable_kernels; 242 for each kernel in append(df_functions,commutator_functions)do 243 solvable_kernels:=delete(kernel,solvable_kernels); 244 if solvable_kernels then 245 246 <<solvable_kernel:= 247 find_solvable_kernel(solvable_kernels, 248 functions_and_constants_list,denominator); 249 if solvable_kernel then 250 <<linear_solve_and_assign(!*ff2a(equation,1),solvable_kernel); 251 depl!*:= 252 delete(assoc(solvable_kernel,depl!*),depl!*); 253 254 successful_message_for(n,"Solved for ",solvable_kernel); 255 goto solved 256 >> 257 else <<not_a_number_message_for(n,"Solving a function", 258 partial_list(solvable_kernels,3)); 259 goto solved>> 260 >> ; 261 262 263 integration_variables:=present_variables; 264 for each kernel in append(linear_functions,commutator_functions)do 265 for each variable in 266 ((if depl_entry then cdr depl_entry) 267 where depl_entry=assoc(kernel,depl!*))do 268 integration_variables:=delete(variable,integration_variables); 269 for each df_function in df_functions do 270 if not (length 271 ((if depl_entry then cdr depl_entry) 272 where depl_entry=assoc(df_function,depl!*))=nr_of_variables) then 273 for each variable in 274 ((if depl_entry then cdr depl_entry) 275 where depl_entry=assoc(df_function,depl!*))do 276 integration_variables:=delete(variable,integration_variables); 277 if try_an_inhomogeneous_integration(n,equation,denominator, 278 df_list,df_terms,integration_variables,nr_of_variables)then goto solved; 279 280 if try_a_differentiation(n,total_used,equation,present_variables, 281 df_terms,linear_functions,commutator_functions) 282 then goto solved; 283 284 write cur_eq_set!*,"(",n,") not solved";terpri!* t; 285 solved: 286 end$ 287 288 289lisp procedure successful_message_for(n,action,kernel); 290 <<write cur_eq_set!*,"(",n,"): ",action; 291 maprin kernel;terpri!*(not !*nat); 292 293 setk(list(cur_eq_set!*,n),0);t>> $ 294 295 296lisp procedure not_a_number_message_for(n,action,kernel); 297 <<write"*** ",cur_eq_set!*,"(",n,"): ",action, 298 " failed:";terpri!* t; 299 write" coefficient not a number for "; 300 maprin kernel;terpri!*(not !*nat); 301 write" Solvable with 'off coefficient_check'"; 302 terpri!* t;t>> $ 303 304 305lisp procedure try_a_homogeneous_integration(n,denominator,df_list); 306 begin scalar solvable_kernel,solvable_kernels,df_kernel; 307 return 308 if null car df_list and 309 (cdr df_list)and length(cdr df_list)=1 310 then 311 if(solvable_kernel:=find_solvable_kernel( 312 solvable_kernels:=list(car car cdr df_list), 313 cdr df_list,denominator))then 314 <<df_kernel:=cadr solvable_kernel; 315 setk(df_kernel,homogeneous_integration_of(solvable_kernel)); 316 depl!*:= 317 delete(assoc(df_kernel,depl!*),depl!*); 318 319 successful_message_for(n,"Homogeneous integration of ",solvable_kernel)>> 320 else not_a_number_message_for(n,"Homogeneous integration", 321 car solvable_kernels) 322 end$ 323 324 325lisp procedure find_solvable_kernel(kernel_list,kc_list,denominator); 326 if !*coefficient_check then 327 first_solvable_kernel(kernel_list,kc_list,denominator) 328 else car kernel_list$ 329 330 331lisp procedure first_solvable_kernel(kernel_list,kc_list,denominator); 332 if kernel_list then 333 (if domainp cdr kc_pair or 334 numberp !*ff2a(cdr kc_pair,denominator) 335 then car kc_pair 336 else first_solvable_kernel(cdr kernel_list,kc_list,denominator)) 337 where kc_pair=assoc(car kernel_list,kc_list)$ 338 339 340lisp procedure homogeneous_integration_of df_term; 341 begin scalar df_function,function_number,dependency_list,integration_list, 342 coefficient_name,bracketname,even_used,odd_used, 343 integration_variable, 344 number_of_integrations,solution,new_dependency_list; 345 346 df_function:=cadr df_term; 347 if not member(car df_function,get(cur_eq_set!*, 'function_list)) 348 or not fixp(function_number:=cadr df_function) 349 or function_number=0 then 350 351 msgpri("PERFORM_HOMOGENEOUS_INTEGRATION: integration of", 352 df_function,"not allowed",nil,t); 353 dependency_list:= 354 ((if depl_entry then cdr depl_entry) 355 where depl_entry=assoc(df_function,depl!*)); 356 if length dependency_list=1 then 357 coefficient_name:=get(cur_eq_set!*, 'constant_operator) 358 else coefficient_name:=car df_function; 359 360 if(bracketname:=get(coefficient_name, 'bracketname))then 361 begin even_used:=get(bracketname, 'even_used); 362 odd_used:=get(bracketname, 'odd_used); 363 end 364 else 365 begin 366 even_used:=get(coefficient_name, 'even_used); 367 odd_used:=get(coefficient_name, 'odd_used); 368 end; 369 integration_list:=cdr cdr df_term; 370 371 if integration_list then integration_variable:=car 372 integration_list else integration_variable:=nil; 373 if integration_variable and(integration_list:=cdr integration_list) 374 and fixp car integration_list then 375 <<number_of_integrations:=car integration_list; 376 integration_list:=cdr integration_list>> 377 else number_of_integrations:=1; 378 if bracketname then 379 380 if function_number>0 then 381 (if even_used+number_of_integrations>get(bracketname, 'even_dimension)then 382 change_dimensions_of(bracketname,even_used+number_of_integrations, 383 get(bracketname, 'odd_dimension))) 384 else 385 (if odd_used+number_of_integrations>get(bracketname, 'odd_dimension)then 386 change_dimensions_of(bracketname,get(bracketname, 'even_dimension), 387 odd_used+number_of_integrations)); 388 389 solution:=nil ./ 1; 390 while integration_variable do 391 begin new_dependency_list:=delete(integration_variable,dependency_list); 392 for i:=0:number_of_integrations-1 do 393 <<solution:=addsq(solution,multsq( 394 if i=0 then 1 ./ 1 else mksq(integration_variable,i), 395 mksq( 396 list(coefficient_name,if function_number>0 then 397 (even_used:=even_used+1)else-(odd_used:=odd_used+1)),1))); 398 if new_dependency_list then 399 depl!*:=(list(coefficient_name,if function_number>0 then even_used 400 else-odd_used) . new_dependency_list) . depl!*; 401 >> ; 402 403 if integration_list then integration_variable:=car 404 integration_list else integration_variable:=nil; 405 if integration_variable and(integration_list:=cdr integration_list) 406 and fixp car integration_list then 407 <<number_of_integrations:=car integration_list; 408 integration_list:=cdr integration_list>> 409 else number_of_integrations:=1 410 411 412 end; 413 solution:=mk!*sq subs2 solution; 414 415 if get(coefficient_name, 'bracketname)then 416 define_used(bracketname,list('list,even_used,odd_used)) 417 else 418 begin 419 put(coefficient_name, 'even_used,even_used); 420 put(coefficient_name, 'odd_used,odd_used); 421 end; 422 return solution 423 end$ 424 425 426lisp procedure split_equation_polynomially(n,total_used,equation, 427 absent_variables); 428 begin scalar polynomial_variables,equations_list; 429 430 polynomial_variables:=absent_variables; 431 if !*polynomial_check then 432 polynomial_variables:=for each variable in polynomial_variables join 433 if polynomialp(equation,variable)then list(variable); 434 435 equations_list:=split_non_linear_form(equation,polynomial_variables); 436 if length equations_list>1 then 437 <<for each pc_pair in cdr equations_list do 438 setk(list(cur_eq_set!*,(total_used:=total_used+1)), 439 mk!*sq((cdr pc_pair) ./ 1)); 440 if car equations_list then 441 setk(list(cur_eq_set!*,(total_used:=total_used+1)), 442 mk!*sq((car equations_list) ./ 1)); 443 write cur_eq_set!*,"(",n,") breaks into ", 444 cur_eq_set!*,"(",get(cur_eq_set!*, 'total_used)+1, 445 "),...,",cur_eq_set!*,"(",total_used,") by "; 446 maprin partial_list(polynomial_variables,5); 447 terpri!*(not !*nat); 448 449 setk(list(cur_eq_set!*,n),0); 450 put(cur_eq_set!*, 'total_used,total_used) 451 >> ; 452 if length equations_list>1 then return t 453 454 455 end$ 456 457 458lisp procedure polynomialp(expression,kernel); 459 if domainp expression then t 460 else((main_variable=kernel or not depends(main_variable,kernel))and 461 polynomialp(lc expression,kernel)and polynomialp(red expression,kernel)) 462 where main_variable=mvar expression$ 463 464 465lisp procedure partial_list(printed_list,nr_of_items); 466 'list . broken_list(printed_list,nr_of_items)$ 467 468lisp procedure broken_list(list,n); 469 if list then if n=0 then '(!.!.!.) 470 else car list . broken_list(cdr list,n-1)$ 471 472 473lisp procedure check_differentiation_sequence(sequence,variable_list); 474 if null sequence then t 475 else if fixp car sequence or 476 member(car sequence,variable_list)then 477 check_differentiation_sequence(cdr sequence,variable_list)$ 478 479 480lisp procedure try_an_inhomogeneous_integration(n,equation,denominator, 481 df_list,df_terms,integration_variables,nr_of_variables); 482 begin scalar solvable_kernel,solvable_kernels,forbidden_functions, 483 df_kernel,inhomogeneous_term; 484 485 for each df_term in df_terms do 486 <<if length 487 ((if depl_entry then cdr depl_entry) 488 where depl_entry=assoc(cadr df_term,depl!*))=nr_of_variables 489 and(check_differentiation_sequence(cdr cdr df_term, 490 integration_variables) 491 or member(cadr df_term,forbidden_functions)) 492 then solvable_kernels:=if member(cadr df_term,forbidden_functions) 493 then list(nil,nil)else df_term . solvable_kernels; 494 forbidden_functions:=(cadr df_term) . forbidden_functions>> ;; 495 496 return 497 if solvable_kernels then 498 if length(solvable_kernels)=1 then 499 if(solvable_kernel:=find_solvable_kernel(solvable_kernels,df_list,denominator)) 500 then 501 if(inhomogeneous_term:=linear_solve(mk!*sq(equation ./ 1),solvable_kernel)) 502 and(not !*polynomial_check or 503 check_polynomial_integration(solvable_kernel,inhomogeneous_term)) 504 then 505 <<df_kernel:=cadr solvable_kernel; 506 setk(df_kernel, 507 inhomogeneous_integration_of(solvable_kernel,inhomogeneous_term)); 508 depl!*:= 509 delete(assoc(df_kernel,depl!*),depl!*); 510 511 successful_message_for(n,"Inhomogeneous integration of ",solvable_kernel)>> 512 else 513 <<write cur_eq_set!*,"(",n,"): Inhomogeneous integration failed: "; 514 terpri!* t; 515 write"inhomogeneous term not polynomial in integration variables"; 516 terpri!* t;t>> 517 else not_a_number_message_for(n,"Inhomogeneous integration", 518 car solvable_kernels) 519 else <<write cur_eq_set!*,"(",n,"): Inhomogeneous integration failed: "; 520 terpri!* t; 521 write"more terms with maximal dependency";terpri!* t;t>> 522 523 524 end$ 525 526lisp procedure check_polynomial_integration(df_term,integration_term); 527 begin scalar numerator,denominator,integration_variables,variable,ok; 528 numerator:=numr simp integration_term; 529 denominator:=denr simp integration_term; 530 integration_variables:= 531 for each argument in cdr cdr df_term join 532 if not fixp argument then list argument; 533 ok:=t; 534 while ok and integration_variables do 535 <<variable:=car integration_variables; 536 ok:=(not depends(denominator,variable)and polynomialp(numerator,variable)); 537 integration_variables:=cdr integration_variables 538 >> ; 539 return ok; 540 end$ 541 542 543lisp procedure inhomogeneous_integration_of(df_term,inhomogeneous_term); 544 begin scalar df_sequence,integration_variables,int_sequence, 545 variable,nr_of_integrations,integration_terms,solution, 546 powers,coefficient,int_factor,solution_term,n,k; 547 df_sequence:=cdr cdr df_term; 548 549 while df_sequence do 550 <<variable:=car df_sequence; 551 df_sequence:=cdr df_sequence; 552 if df_sequence and fixp car df_sequence then 553 <<nr_of_integrations:=car df_sequence; 554 df_sequence:=cdr df_sequence>> 555 else nr_of_integrations:=1; 556 integration_variables:=variable . integration_variables; 557 int_sequence:=(variable . nr_of_integrations) . int_sequence 558 >> ; 559 integration_terms:=split_non_linear_form(numr simp inhomogeneous_term, 560 integration_variables); 561 integration_terms:=(nil . car integration_terms) . 562 cdr integration_terms; 563 564 565 solution:=nil ./ 1; 566 for each term in integration_terms do 567 <<powers:=car term;coefficient:=cdr term; 568 int_factor:=1;solution_term:=1 ./ 1; 569 for each integration in int_sequence do 570 <<variable:=car integration;k:=cdr integration; 571 n:=(if power then cdr power else 0)where power=assoc(variable,powers); 572 573 for i:=1:k do int_factor:=(n+i)*int_factor; 574 solution_term:=multsq(solution_term,mksq(variable,n+k)) 575 >> ; 576 solution_term:=multsq(solution_term,coefficient ./ int_factor); 577 solution:=addsq(solution,solution_term) 578 >> ; 579 solution:=multsq(solution,1 ./ denr simp inhomogeneous_term); 580 solution:=mk!*sq subs2 addsq(solution,simp homogeneous_integration_of df_term); 581 return solution 582 end$ 583 584 585lisp procedure try_a_differentiation(n,total_used,equation,present_variables, 586 df_terms,linear_functions,commutator_functions); 587 begin scalar differentiations_list,polynomial_order; 588 589 590 present_variables:=for each variable in present_variables collect 591 (variable . nil . 0); 592 593 for each kernel in df_terms do 594 for each variable in 595 ((if depl_entry then cdr depl_entry) 596 where depl_entry=assoc(cadr(kernel),depl!*))do 597 598 rplacd(entry,kernel . (cddr entry+1)) 599 where entry=assoc(variable,present_variables);; 600 601 for each kernel in linear_functions do 602 for each variable in 603 ((if depl_entry then cdr depl_entry) 604 where depl_entry=assoc(kernel,depl!*))do 605 606 rplacd(entry,kernel . (cddr entry+1)) 607 where entry=assoc(variable,present_variables);; 608 609 for each kernel in commutator_functions do 610 for each variable in 611 ((if depl_entry then cdr depl_entry) 612 where depl_entry=assoc(kernel,depl!*))do 613 614 rplacd(entry,nil . (cddr entry+1)) 615 where entry=assoc(variable,present_variables);; 616 617 differentiations_list:= 618 for each entry in present_variables join 619 if cadr entry and cddr entry=1 and 620 (polynomial_order:=get_polynomial_order( 621 linear_solve(mk!*sq(equation ./ 1),cadr entry),car entry)) 622 then list(car entry . cadr entry . (polynomial_order+1)); 623 return 624 if differentiations_list then 625 if !*allow_differentiation then 626 <<for each entry in differentiations_list do 627 setk(list(cur_eq_set!*,(total_used:=total_used+1)), 628 mk!*sq simpdf list(mk!*sq(equation ./ 1), 629 car entry,cddr entry)); 630 write cur_eq_set!*,"(",n,"): Generation of ", 631 cur_eq_set!*,"(",get(cur_eq_set!*, 'total_used)+1, 632 "),...,",cur_eq_set!*,"(",total_used, 633 ") by differentiation w.r.t. "; 634 terpri!* t; 635 maprin partial_list(for each entry in differentiations_list collect 636 list('list,car entry,cddr entry),10); 637 terpri!*(not !*nat); 638 put(cur_eq_set!*, 'total_used,total_used); 639 t 640 >> 641 else << 642 write"*** ",cur_eq_set!*,"(",n, 643 "): Generation of new equations by differentiation possible."; 644 terpri!* t;write" Solvable with 'on allow_differentiation'"; 645 terpri!* t;t>> 646 647 648 end$ 649 650 651lisp procedure get_polynomial_order(expression,variable); 652 if not depends(denr(expression:=simp expression),variable)and 653 (not !*polynomial_check or polynomialp(numr expression,variable))then 654 begin scalar kord!*; 655 setkorder list !*a2k variable; 656 expression:=reorder numr expression; 657 return if mvar expression=variable then ldeg expression else 0; 658 end$ 659 660 661algebraic procedure integrate_equations(m,n); 662 for i:=m:n do integrate_equation(i)$ 663 664 665lisp operator integrate_exceptional_equation; 666lisp procedure integrate_exceptional_equation(n); 667 integrate_equation(n) 668 where 669 !*coefficient_check=nil, 670 !*polynomial_check=nil, 671 !*allow_differentiation=t$ 672 673 674 675lisp operator auto_solve; 676lisp procedure auto_solve nr_list; 677 begin scalar total,old_total,to_do,unsolved,old_unsolved,stuck; 678 total:=old_total:=get(cur_eq_set!*, 'total_used); 679 to_do:=if fixp nr_list then list nr_list 680 else if car nr_list= 'list then cdr nr_list 681 else nr_list; 682 while not stuck and to_do do begin 683 for each eq_nr in to_do do 684 <<integrate_equation eq_nr; 685 if cadr assoc(list(cur_eq_set!*,eq_nr),get(cur_eq_set!*, 'kvalue))neq 0 686 then unsolved:=eq_nr . unsolved>> ; 687 total:=get(cur_eq_set!*, 'total_used); 688 if total=old_total and unsolved and unsolved=old_unsolved then stuck:=t 689 else <<old_unsolved:=unsolved; 690 to_do:=reverse unsolved; 691 unsolved:=nil; 692 to_do:=append(for eq_nr:=old_total+1:total collect eq_nr,to_do); 693 old_total:=total>> 694 end; 695 if stuck then return 'list . reverse unsolved 696 else <<terpri();write"Successful integration of all equations";terpri()>> ; 697 end$ 698 699lisp operator show_equation; 700lisp procedure show_equation n; 701 begin scalar equation,total_used,function_list; 702 if null(total_used:=get(cur_eq_set!*, 'total_used))or 703 n>total_used then 704 705 msgpri("SHOW_EQUATION: properly initialize", 706 cur_eq_set!*,nil,nil,t); 707 if(equation:=assoc(list(cur_eq_set!*,n), 708 get(cur_eq_set!*, 'kvalue)))then 709 begin 710 equation:=setk(list(cur_eq_set!*,n),aeval cadr equation); 711% Modified to keep into account the change from varpri to assignpri 712% in Reduce 3.5 713% varpri(equation,list('setk,mkquote list(cur_eq_set!*,n), 714% mkquote equation), 'only); 715 assgnpri(equation,list(list(cur_eq_set!*,n)), 'only); 716 function_list:=get_recursive_kernels(numr simp equation, 717 get(cur_eq_set!*, 'function_list)); 718 if function_list then 719 <<terpri!* t;write"Functions occurring:";terpri!* t; 720 for each fn in function_list do 721 <<maprin(fn . 722 ((if depl_entry then cdr depl_entry) 723 where depl_entry=assoc(fn,depl!*)));terpri!*(not !*nat)>> 724 >> 725 else terpri!* nil 726 end 727 end$ 728 729 730algebraic procedure show_equations(m,n); 731 for i:=m:n do show_equation i$ 732 733 734lisp operator functions_used,put_functions_used, 735 equations_used,put_equations_used; 736 737 738lisp procedure functions_used function_name; 739 list('list,get(function_name, 'even_used),get(function_name, 'odd_used))$ 740 741 742lisp procedure put_functions_used(function_name,even_used,odd_used); 743 begin 744 if not fixp even_used or even_used<0 or 745 not fixp odd_used or odd_used<0 then 746 747 msgpri("PUT_FUNCTIONS_USED: used functions number invalid", 748 nil,nil,nil,t); 749 put(function_name, 'even_used,even_used); 750 put(function_name, 'odd_used,odd_used); 751 end$ 752 753 754lisp procedure equations_used; 755 get(cur_eq_set!*, 'total_used)$ 756 757 758lisp procedure put_equations_used(n); 759 if not fixp n or n<0 then 760 761 msgpri("PUT_EQUATIONS_USED: used equation number invalid",nil,nil,nil,t) 762 else put(cur_eq_set!*, 'total_used,n)$ 763 764 765lisp operator df_acts_as_derivation_on; 766 767lisp procedure df_acts_as_derivation_on operator_name; 768 begin 769 put(operator_name, 'dfform, 'df_as_derivation); 770 end$ 771 772 773lisp procedure df_as_derivation(kernel,variable,power); 774 begin scalar left_part,right_part,argument,derivative; 775 if power neq 1 then 776 777 msgpri("DF_AS_DERIVATION:",kernel,"must occur linearly",nil,t); 778 left_part:=list car kernel; 779 right_part:=cdr kernel; 780 derivative:=nil . 1; 781 while right_part do 782 <<argument:=car right_part; 783 right_part:=cdr right_part; 784 derivative:=addsq(derivative, 785 simp append(reverse left_part, 786 list('df,argument,variable) . right_part)); 787 left_part:=argument . left_part; 788 >> ; 789 return derivative; 790 end$ 791 792 793lisp operator listlength$ 794 795lisp procedure listlength l; 796 listpri_depth!*:=l$ 797 798 799symbolic procedure listpri l; 800 begin 801 scalar orig,split,u; 802 u:=l; 803 l:=cdr l; 804 prin2!* get('!*lcbkt!*, 'prtch); 805 orig:=orig!*; 806 orig!*:=if posn!*<18 then posn!* else orig!*+3; 807 if null l then go to b; 808 split:=treesizep(l,listpri_depth!*); 809 a: maprint(negnumberchk car l,0); 810 l:=cdr l; 811 if null l then go to b; 812 oprin '!*comma!*; 813 if split then terpri!* t; 814 go to a; 815 b: prin2!* get('!*rcbkt!*, 'prtch); 816 orig!*:=orig; 817 return u 818 end$ 819 820endmodule; 821end; 822 823 824