1module desir; % Special case differential equation solver. 2 3% Redistribution and use in source and binary forms, with or without 4% modification, are permitted provided that the following conditions are met: 5% 6% * Redistributions of source code must retain the relevant copyright 7% notice, this list of conditions and the following disclaimer. 8% * Redistributions in binary form must reproduce the above copyright 9% notice, this list of conditions and the following disclaimer in the 10% documentation and/or other materials provided with the distribution. 11% 12% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 13% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 14% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 15% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 16% CONTRIBUTORS 17% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 18% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 19% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 20% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 21% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 22% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 23% POSSIBILITY OF SUCH DAMAGE. 24% 25 26 27create!-package('(desir),'(solve)); 28 29% *************************************************************** 30% * * 31% * DESIR * 32% * ===== * 33% * * 34% * SOLUTIONS FORMELLES D'EQUATIONS DIFFERENTIELLES * 35% * * 36% * LINEAIRES ET HOMOGENES * 37% * * 38% * AU VOISINAGE DE POINTS SINGULIERS REGULIERS ET IRREGULIERS * 39% * * 40% *************************************************************** 41% 42% Differential linear homogenous Equation Solutions in the 43% neighbourhood of Irregular and Regular points 44% 45% Version 3.3 - Novembre 1993 46% 47% 48% Groupe de Calcul Formel de Grenoble 49% laboratoire LMC 50% 51% E-mail: dicresc@imag.fr 52 53 54 55 56 57 58% This software enables the basis of formal solutions to be computed 59% for an ordinary homogeneous differential equation with polynomial 60% coefficients over Q of any order, in the neighbourhood of zero 61% (regular or irregular singular point, or ordinary point ). 62% Tools have been added to deal with equations with a polynomial 63% right-hand side, parameters and a singular point not to be found at 64% zero. 65% 66% This software can be used in two ways : * direct ( DELIRE procedure) 67% * interactive ( DESIR procedure) 68% 69% The basic procedure is the DELIRE procedure which enables the 70% solutions of a linear homogeneous differential equation to be 71% computed in the neighbourhood of zero. 72% 73% The DESIR procedure is a procedure without argument whereby DELIRE 74% can be called without preliminary treatment to the data, that is to 75% say, in an interactive autonomous way. This procedure also proposes 76% some transformations on the initial equation. This allows one to 77% start comfortably with an equation which has a non zero singular 78% point, a polynomial right-hand side and parameters. 79% 80% This document is a succint user manual. For more details on the 81% underlying mathematics and the algorithms used, the reader can refer 82% to : 83% 84% E. Tournier : Solutions formelles d'equations differentielles - 85% Le logiciel de calcul formel DESIR. 86% These d'Etat de l'Universite Joseph Fourier (Grenoble, Avr. 87). 87% 88% He will find more precision on use of parameters in : 89% 90% F. Richard-Jung : Representation graphique de solutions 91% d'equations differentielles dans le champ complexe. 92% These de l'Universite Louis Pasteur (Strasbourg - septembre 88). 93 94 95 96% 97 98 99% ************************** 100% * * 101% * FORMS OF SOLUTIONS * 102% * * 103% ************************** 104 105 106 107 108% We have tried to represent solutions in the simplest form possible. 109% For that, we have had to choose different forms according to the 110% complexity of the equation (parameters) and the later use we shall 111% have of these solutions. 112% 113% "general solution" = {......, { split_sol , cond },....} 114% ------------------ 115% 116% cond = list of conditions or empty list (if there is no condition) 117% that parameters have to verify such that split_sol is in 118% the basis of solutions. In fact, if there are parameters, 119% basis of solutions can have different expressions according 120% to the values of parameters. ( Note : if cond={}, the list 121% "general solution" has one element only. 122% 123% split_sol = { q , ram , polysol , r } 124% ( " split solution " enables precise information on the 125% solution to be obtained immediately ) 126% 127% The variable in the differential operator being x, solutions are 128% expressed with respect to a new variable xt, which is a fractional 129% power of x, in the following way : 130% 131% q : polynomial in 1/xt with complex coefficients 132% ram : xt = x**ram (1/ram is an integer) 133% polysol : polynomial in log(xt) with formal series in xt 134% coefficients 135% r : root of a complex coefficient polynomial ("indicial 136% equation"). 137% 138% 139% qx r*ram 140% "standard solution" = e x polysolx 141% ----------------- 142% qx and polysolx are q and polysol expressions in which xt has 143% been replaced by x**ram 144% 145% N.B. : the form of these solutions is simplified according to the 146% nature of the point zero. 147% - if 0 is a regular singular point : the series appearing in polysol 148% are convergent, ram = 1 and q = 0. 149% - if 0 is a regular point, we also have : polysol is constant in 150% log(xt) (no logarithmic terms). 151 152 153% 154 155 156% *********************************** 157% * * 158% * INTERACTIVE USE * 159% * * 160% *********************************** 161% 162 163% Modification of the "deg" function. 164 165symbolic procedure deg(u,kern); 166 <<u := simp!* u; tstpolyarg(denr u,kern); numrdeg(numr u,kern)>> 167 where dmode!* = gdmode!*; 168 169%symbolic procedure deg(u,kern); 170% begin scalar x,y; 171% u := simp!* u; 172% y := denr u; 173% tstpolyarg(y,u); 174% u := numr u; 175% kern := !*a2k kern; 176% if domainp u then return 0 177% else if mvar u eq kern then return !*f2a ldeg u; 178% x := setkorder list kern; 179% u := reorder u; 180%% if not(mvar u eq kern) then u := nil else u := ldeg u; 181% if not(mvar u eq kern) then u := 0 else u := ldeg u; 182% setkorder x; 183% return !*f2a u 184% end; 185 186fluid '(!*precise !*trdesir); 187 188switch trdesir; 189 190global '(multiplicities!*); 191 192flag('(multiplicities!*),'share);% Since SOLVE not loaded when file 193 % compiled. 194 195!*precise := nil; % Until we understand the interactions. 196 197algebraic; 198 199procedure desir ; 200%===============; 201% 202% Calcul des solutions formelles d'une equation differentielle lineaire 203% homogene de maniere interactive. 204% La variable dans cette equation est obligatoirement x. 205% ----------------- 206% x et z doivent etre des variables atomiques. 207% 208% La procedure demande l'ordre et les coefficients de l'equation, le 209% nom des parametres s'il y en a, puis si l'on souhaite une 210% transformation de cette equation et laquelle ( par exemple, ramener 211% un point singulier a l'origine - voir les procedures changehom, 212% changevar, changefonc - ). 213% 214% Cette procedure ECRIT les solutions et RETOURNE une liste de terme 215% general { lcoeff, {....,{ solution_generale },....}}. Le nombre 216% d'elements de cette liste est lie au nombre de transformations 217% demandees : 218% * lcoeff : liste des coefficients de l'equation differentielle 219% * solution_generale : solution ecrite sous la forme generale 220begin 221 scalar k,grille,repetition,lcoeff,param,n,ns,solutions,lsol ; 222 integer j; 223 if (repetition neq non ) and (repetition neq nonon ) then 224 << write 225 " ATTENTION : chaque donnee doit etre suivie de ; ou de $" ; 226 repetition:=nonon ; 227 >> ; 228 solutions:={}; 229 lsol:={} ; 230 % lecture des donnees ; 231 lcoeff:= lectabcoef(); 232 param:=second lcoeff; 233 lcoeff:=first lcoeff; 234 continue:=oui; 235 write "transformation ? (oui;/non;) "; 236 ok:=xread(nil); 237 while continue eq oui do 238 << 239 if ok eq oui then <<lcoeff:=transformation(lcoeff,param); 240 param:=second lcoeff; 241 lcoeff:=first lcoeff; 242 >>; 243 244 write "nombre de termes desires pour la solution ?" ; 245 k:=xread(nil) ; 246 if k neq 0 then 247 << 248 grille:=1 ; 249 if repetition neq non and lisp !*trdesir then 250 << write " "; 251 write "a chaque etape le polygone NRM sera visualise par la ", 252 "donnee des aretes modifiees , sous la forme :" ; 253 write " " ; 254 write 255 " ARETE No i : coordonnees de l'origine gauche, pente,", 256 " longueur " ; 257 >> ; 258 write " " ; 259 on div ; 260 261 on gcd ; 262 solutions:= delire(x,k,grille,lcoeff,param); 263 ns:=length solutions; n:=length lcoeff -1; 264 if ns neq 0 then 265 << write "LES ",ns," SOLUTIONS CALCULEES SONT LES SUIVANTES"; 266 j:=1; 267 for each elt in solutions do 268 << write " " ; write " ==============" ; 269 write " SOLUTION No ",j ; 270 write " ==============" ; 271 sorsol(elt); 272 j:=j+1; 273 >> ; 274 >>; 275 off div ; 276 if ns neq n then write n-ns," solutions n'ont pu etre calculees"; 277 repetition:=non ; 278 lsol:= append(lsol,{{lcoeff,solutions}}) ; 279 write "voulez-vous continuer ? "; 280 write 281 "'non;' : la liste des solutions calculees est affichee (sous"; 282 write " forme generalisee)."; 283 write "'non$' : cette liste n'est pas affichee."; 284 continue:=xread(nil); ok:=oui; 285 >> 286 else 287 continue:=non; 288 >>; 289 290 return lsol ; 291end; 292 293procedure solvalide(solutions,solk,k) ; 294%==================================== ; 295% 296% Verification de la validite de la solution numero solk dans la liste 297% solutions : {lcoeff,{....,{solution_generale},....}}. 298% On reporte la solution dans l'equation : le resultat a en facteur un 299% polynome en xt qui doit etre de degre > une valeur calculee en 300% fonction de k, nombre de termes demandes dans le developpement 301% asymptotique. Ne peut etre utilisee si la solution numero solk est 302% liee a une condition. 303% 304% ECRIT et RETOURNE l'evaluation de ce report. 305begin 306 scalar z,lcoeff,sol,essai,qx,gri,r,coeff1,d,zz; 307 integer j; 308 lcoeff:=first solutions; 309 sol:=part(second solutions,solk); 310 if length sol > 1 311 then write("presence de solutions conditionnelles :", 312 " cette procedure ne peut pas etre appelee.") 313 else 314 << z:=first sol; 315 essai:=first z; qx:=first essai; 316 essai:=rest essai; 317 gri:= first essai; 318 sol:=second essai; r:=third essai; 319 essai:=second z ;if length(essai)>0 then 320 write "presence d'une condition : cette procedure ne peut pas etre 321 appelee." 322 else 323 <<%calcul de la valuation theorique du polynome reste 324 coeff1:=for each elt in lcoeff collect 325 sub(x=xt**(1/gri),elt); 326 if qx neq 0 then <<d:=solvireg(coeff1,qx,xt); 327 coeff1:=changefonc(coeff1,xt,!&phi,e**qx*!&phi); 328 >>; 329 d:=altmin(coeff1,xt)-d; 330 331 qx:=sub(xt=x**gri,qx); 332 sol:=sub(lambd=r,sol); 333 sol:=e**qx*x**(r*gri)*sub(xt=x**gri,sol); 334 write "La solution numero ",solk," est ",sol; 335 write "La partie reguliere du reste est de l'ordre de x**(", 336 gri*(k+1+d+r),")"; 337 z:=0; 338 for each elt in lcoeff do 339 << z:=z+elt*df(sol,x,j);j:=j+1;>>; 340 zz:=e**(-qx)*x**(-r*gri)*z; 341 zz:=sub(x=xt**(1/gri),zz); 342 on rational; 343 write "Si on reporte cette solution dans l'equation, le terme ", 344 "significatif du reste"," est : ", 345 e**qx*x**(r*gri)*sub(xt=x**gri,valterm(zz,xt)); 346 off rational; 347 return z ; 348 >> ; 349 >>; 350end; 351 352procedure solvireg(lcoeff,q,x); 353%=============================; 354begin scalar f; 355 integer j,n; 356 depend !&y,x; 357 depend !&phi,x; 358 l:=lcoeff; 359 while l neq {} do 360 <<f:=f+(first l)*df(!&y,x,j);j:=j+1;l:=rest l>>; 361 n:=length(lcoeff); 362 let !&y=e**q*!φ 363 for j:=1:n do f:=sub(df(!&phi,x,j)=zz**j,f); 364 f:=sub(!&phi=1,f); 365 clear !&y; 366 nodepend !&y,x; 367 nodepend !&phi,x; 368 return deg(den(f),x); 369end; 370 371procedure altmin(lcoeff,x); 372%=========================; 373begin 374integer j,min,d; 375min:=deg(valterm(first lcoeff,x),x); 376for each elt in rest lcoeff do << 377 j:=j+1; 378 d:=deg(valterm(elt,x),x); 379 if d-j<min then min:=d-j;>>; 380return min; 381end; 382 383procedure valterm(poly,x); 384%=========================; 385%retourne le terme de plus bas degre de poly; 386begin 387scalar l,elt; 388integer j; 389l:=coeff(poly,x); 390elt:=first l; 391while (elt=0) and (rest(l) neq {}) do 392 <<j:=j+1;l:=rest l; elt:=first l>>; 393return elt*x**j; 394end; 395 396procedure standsol(solutions) ; 397%============================== ; 398% 399% PERMET d'avoir l'expression simplifiee de chaque solution a partir de 400% la liste des solutions {lcoeff,{....,{solution_generale},....}}, qui 401% est retournee par DELIRE ou qui est un des elements de la liste 402% retournee par DESIR. 403% 404% RETOURNE une liste de 3 elements : { lcoeff , solstand, solcond } 405% * lcoef = liste des coefficients de l'equation differentielle 406% * solstand = liste des solutions sous la forme standard 407% * solcond = liste des solutions conditionnelles n'ayant pu etre 408% mises sous la forme standard. Ces solutions restent 409% sous la forme generales 410% 411% Cette procedure n'a pas de sens pour les solutions "conditionnelles". 412% Pour celles-ci, il est indispensable de donner une valeur aux 413% parametres, ce que l'on peut faire, soit en appelant la procedure 414% SORPARAM, qui ecrit et retourne ces solutions dans la forme standard, 415% soit en appelant la procedure SOLPARAM qui les retourne dans la forme 416% generale. 417begin 418 scalar z,lcoeff,sol,solnew,solcond,essai,qx,gri,r; 419 integer j; 420 solnew:={}; 421 solcond:= {} ; 422 lcoeff:=first solutions; 423 for each elt in second solutions do 424 if length elt > 1 then solcond:=append(solcond,{elt}) 425 else 426 << z:=first elt; 427 essai:=first z; 428 qx:=first essai; 429 essai:=rest essai; 430 gri:= first essai; 431 qx:=sub(xt=x**gri,qx); 432 sol:=second essai; r:=third essai; 433 essai:=second z ; 434 if length(essai)>0 then solcond:=append(solcond,{elt}) 435 else 436 << sol:=sub(lambd=r,sol); 437 sol:=e**qx*x**(r*gri)*sub(xt=x**gri,sol); 438 solnew:=append(solnew,{sol}); 439 >> ; 440 >>; 441return {lcoeff,solnew,solcond}; 442end; 443 444procedure sorsol(sol); 445%===================== 446% 447% ecriture, sous forme standard, de la solution sol donnee sous la forme 448% generale, avec enumeration des differentes conditions (s'il y a lieu). 449% 450begin 451 scalar essai,qx,gri,sol,r; 452 nonnul:=" non nul"; 453 entnul:=" nul"; 454 nonent:=" non entier" ; 455 entpos:= " entier positif" ; 456 entneg:= " entier negatif" ; 457 458 for each z in sol do 459 << essai:=first z; 460 qx:=first essai; 461 essai:=rest essai; 462 gri:= first essai; 463 qx:=sub(xt=x**gri,qx); 464 sol:=second essai; 465 r:=third essai; 466 essai:=second z ; 467 sol:=sub(xt=x**gri,sol); 468 if length(essai)>0 then 469 <<if deg(num sol,lambd)=0 then 470 <<write e**qx*x**(r*gri)*sol ; 471 write "Si : "; 472 for each w in essai do 473 if (length(w)=2 or not lisp !*trdesir) then 474 write first w,second w 475 else 476 << write (first w,second w,third w); 477 w:=rest rest rest w; 478 for each w1 in w do 479 write (" +-",w1); 480 >> 481 >> 482 >> 483 else 484 << sol:=sub(lambd=r,sol); 485 write e**qx*x**(r*gri)*sol; 486 >>; 487 >>; 488clear nonnul,entnul,nonent,entpos,entneg; 489end; 490 491procedure changehom(lcoeff,x,secmembre,id); 492%======================================== 493% 494% derivation d'une equation avec second membre. 495% * lcoeff : liste des coefficients de l'equation 496% * x : variable 497% * secmembre : second membre 498% * id : ordre de la derivation 499% 500% retourne la liste des coefficients de l'equation derivee 501% permet de transforme une equation avec second membre polynomial en une 502% equation homogene en derivant id fois, id = degre(secmembre) + 1. 503% 504begin scalar l,fct,cf,n; 505 integer j; 506 depend !&y,x; 507 fct:=secmembre; 508 l:=lcoeff; 509 while l neq {} do 510 <<fct:=fct+(first l)*df(!&y,x,j);j:=j+1;l:=rest l>>; 511 fct:=df(fct,x,id); 512 n:=length(lcoeff)+id; 513 for j:=1:n do fct:=sub(df(!&y,x,j)=zz**j,fct); 514 fct:=sub(!&y=1,fct); 515 cf:=coeff(fct,zz); 516 j:=0; 517 if lisp !*trdesir then 518 for each elt in cf do <<write "a(",j,") = ", elt;j:=j+1;>>; 519 nodepend !&y,x; 520 return cf; 521end; 522 523procedure changevar(lcoeff,x,v,fct); 524%================================= 525% 526% changement de variable dans l'equation homogene definie par la liste, 527% lcoeff, de ses coefficients : 528% l'ancienne variable x et la nouvelle variable v sont liees par la 529% relation x = fct(v) 530% 531% retourne la liste des coefficients en la variable v de la nouvelle 532% equation 533% Exemples d'utilisation : 534% - translation permettant de ramener une singularite rationnelle a 535% l'origine 536% - x = 1/v ramene l'infini en 0. 537% 538begin scalar f,cf; 539 integer j,n; 540 depend !&y,x; 541 l:=lcoeff; 542 while l neq {} do 543 <<f:=f+(first l)*df(!&y,x,j);j:=j+1;l:=rest l>>; 544 n:=length(lcoeff); 545 f:=change(!&y,x,v,fct,f,n); 546 for j:=1:n do f:=sub(df(!&y,v,j)=zz**j,f); 547 f:=sub(!&y=1,f); 548 cf:=coeff(num(f),zz); 549 j:=0; 550 if lisp !*trdesir then 551 for each elt in cf do <<write "a(",j,") = ", elt;j:=j+1;>>; 552 nodepend !&y,x; 553 return cf; 554 end; 555 556procedure changefonc(lcoeff,x,q,fct); 557%================================ 558% 559% changement de fonction inconnue dans l'equation homogene definie par 560% la liste lcoeff de ses coefficients : 561% * lcoeff : liste des coefficients de l'equation initiale 562% * x : variable 563% * q : nouvelle fonction inconnue 564% * fct : y etant la fonction inconnue y = fct(q) 565% 566% retourne la liste des coefficients de la nouvelle equation 567% 568% Exemple d'utilisation : permet de calculer, au voisinage d'une 569% singularite irreguliere, l'equation reduite associee a l'une des 570% pentes (polygone de Newton ayant une pente nulle de longueur non 571% nulle). Cette equation fournit de nombreux renseignements sur la 572% serie divergente associee. 573% 574begin scalar f,cf; 575 integer j,n; 576 depend !&y,x; 577 depend q,x; 578 l:=lcoeff; 579 while l neq {} do 580 <<f:=f+(first l)*df(!&y,x,j);j:=j+1;l:=rest l>>; 581 n:=length(lcoeff); 582 let !&y=fct; 583 for j:=1:n do f:=sub(df(q,x,j)=zz**j,f); 584 f:=sub(q=1,f); 585 cf:=coeff(num(f),zz); j:=1; 586 if lisp !*trdesir then 587 for each elt in cf do <<write "a(",j,") = ", elt;j:=j+1;>>; 588 clear !&y; 589 nodepend !&y,x; 590 nodepend q,x; 591 return cf; 592end; 593 594procedure sorparam(solutions,param); 595%================================== 596% 597% procedure interactive d'ecriture des solutions evaluees : la valeur 598% des parametres est demandee. 599% * solutions : {lcoeff,{....,{solution_generale},....}} 600% * param : liste des parametres; 601% 602% retourne la liste formee des 2 elements : 603% * liste des coefficients evalues de l'equation 604% * liste des solutions standards evaluees pour les valeurs des 605% parametres 606% 607begin 608 scalar essai,sec,qx,gri,sol,sol1,sol2,r,solnew,coefnew,omega,omegac; 609 integer j,iparam; 610 solnew:={}; 611 iparam:=length param; 612 if iparam=0 613 then rederr "La liste des parametres est vide : utiliser STANDSOL"; 614 array parm(iparam),parmval(iparam); 615 j:=1; 616 for each elt in param do 617 << write "donner la valeur du parametre ", elt; 618 parm(j):=elt;parmval(j):=xread(nil);j:=j+1; 619 >>; 620 j:=1; 621 for each elt in second solutions do 622 << for each z in elt do 623 << essai:=first z; 624 qx:=first essai; 625 essai:=rest essai; 626 gri:= first essai; 627 qx:=sub(xt=x**gri,qx); 628 sol1:=second essai; 629 r:=third essai; 630 essai:=second z ; 631 if essai ={} then 632 << sol:=e**qx*x**(r*gri)*sub(xt=x**gri,sol1); 633 for j:=1:iparam do sol:=sub(parm(j)=parmval(j),sol); 634 >> 635 else <<sol2:=sorparamcond(essai,iparam,qx,gri,r,sol1); 636 if sol2 neq 0 then sol:=sol2; >>; 637 >>; 638 write " " ; write " ==============" ; 639 write " SOLUTION No ",j ; 640 write " ==============" ; 641 if sol neq 0 then 642 <<write sol; solnew:=append(solnew,{sol})>> 643 else write "solution non calculee"; 644 j:=j+1; 645 >> ; 646 coefnew:= for each elt in first solutions collect 647 begin scalar cof; 648 cof:=elt ; 649 for j:=1:iparam do cof:=sub(parm(j)=parmval(j),cof); 650 return cof 651 end; 652 clear parm,parmval; 653 return { coefnew, solnew }; 654end; 655 656procedure sorparamcond(essai,iparam,qx,gri,r,sol1); 657%=================================================; 658begin 659scalar sol,sec,omega,omegac; 660 essai:=first essai ; 661 omega:=first essai; 662 sec:= second essai ; 663 for j:=1:iparam do omega:=sub(parm(j)=parmval(j),omega); 664 omegac:=append(coeff(omega,i),{0}); 665 on rounded; 666 if not numberp(first omegac) or not numberp(second omegac) 667 then rederr list("Les valeurs donnees aux parametres ne", 668 "permettent pas de choisir parmi les solutions conditionnelles."); 669 off rounded; 670 % il ne faut traiter qu'une seule fois la solution 671 if sec=nonnul then 672 if omega neq 0 then 673 << sol:=e**qx*x**(r*gri)*sub(xt=x**gri,sol1); 674 for j:=1:iparam do sol:=sub(parm(j)=parmval(j),sol); 675 >>; 676 if sec= entnul then 677 if omega=0 then 678 << sol:=e**qx*x**(r*gri)*sub(xt=x**gri,sol1); 679 for j:=1:iparam do sol:=sub(parm(j)=parmval(j),sol); 680 >>; 681 if sec=nonent then 682 if not fixp(omega) then 683 << sol:=e**qx*x**(r*gri)*sub(xt=x**gri,sol1); 684 for j:=1:iparam do sol:=sub(parm(j)=parmval(j),sol); 685 >>; 686 if sec=entpos then 687 if fixp(omega) and (omega>0) then 688 << sol:=e**qx*x**(r*gri)*sub(xt=x**gri,sol1); 689 for j:=1:iparam do sol:=sub(parm(j)=parmval(j),sol); 690 >>; 691 692 if sec=entneg then 693 if fixp(omega) and (omega<0) then 694 << sol:=e**qx*x**(r*gri)*sub(xt=x**gri,sol1); 695 for j:=1:iparam do sol:=sub(parm(j)=parmval(j),sol); 696 >>; 697 698 if deg(num sol,lambd) neq 0 then 699 << sol:=sub(lambd=r,sol); 700 for j:=1:iparam do sol:=sub(parm(j)=parmval(j),sol); 701 >>; 702return sol; 703end; 704 705procedure solparam(solutions,param,valparam); 706%=========================================== 707% 708% Cette procedure evalue, pour les valeurs des parametres donnees dans 709% valparam les solutions generalisees et les retourne sous forme 710% generalisee. 711% 712% * solutions : {lcoeff,{....,{solution_generale},....}} 713% * param : liste des parametres; 714% * valparam : liste des valeurs des parametres 715% 716% retourne la liste formee des 2 elements : 717% * liste des coefficients evalues de l'equation 718% * liste des solutions sous la forme generalisee evaluees pour les 719% valeurs des parametres 720% 721begin 722 scalar essai,sol,sol1,solg,solnew, coefnew; 723 integer j,iparam; 724 solnew:={}; 725 iparam:=length param; 726if iparam=0 727 then rederr "La liste des parametres est vide : utiliser STANDSOL"; 728 array parm(iparam),parmval(iparam); 729 j:=1; 730 for each elt in param do 731 << parm(j):=elt ; j:=j+1 >>; 732 j:=1; 733 for each elt in valparam do 734 << parmval(j):=elt ; j:=j+1 >>; 735 for each elt in second solutions do 736 << for each z in elt do 737 << solg:=first z; 738 essai:=second z ; 739 if essai ={} then 740 sol1:=solg 741 else sol1:=solparamcond(essai,iparam,solg); 742 if sol1 neq {} then 743 << essai:=rest(rest(sol1)) ; r:=second essai; 744 if deg(num(sol:=first(essai)),lambd) neq 0 then 745 << sol:=sub(lambd=r,sol); 746 for j:=1:iparam do 747 sol:=sub(parm(j)=parmval(j),sol); 748 >>; 749 sol1:={first(sol1), second(sol1),sol,r}; 750 solnew:=append(solnew,{{{sol1,{}}}}); 751 >> ; 752 >>; 753 >> ; 754 coefnew:= for each elt in first solutions collect 755 begin scalar cof; 756 cof:=elt ; 757 for j:=1:iparam do cof:=sub(parm(j)=parmval(j),cof); 758 return cof 759 end; 760 clear parm,parmval; 761 return { coefnew, solnew }; 762end; 763 764 765procedure solparamcond(essai,iparam,solg); 766%========================================; 767begin 768scalar sec,sol1,sol,omega,omegac; 769 essai:=first essai ; 770 omega:=first essai; 771 sec:= second essai ; 772 for j:=1:iparam do omega:=sub(parm(j)=parmval(j),omega); 773 omegac:=append(coeff(omega,i),{0}); 774 on rounded; 775 if not numberp(first omegac) or not numberp(second omegac) 776 then rederr list("Les valeurs donnees aux parametres", 777 "ne permettent pas de choisir parmi les solutions conditionnelles."); 778 off rounded; 779 % il ne faut traiter qu'une seule fois la solution 780 sol1:={}; 781 if sec= nonnul then 782 if omega neq 0 then 783 sol1:= for each elem in solg collect 784 begin 785 sol:=elem ; 786 for j:=1:iparam do 787 sol:=sub(parm(j)=parmval(j),sol); 788 return sol 789 end ; 790 if sec= entnul then 791 if omega=0 then 792 sol1:= for each elem in solg collect 793 begin 794 sol:=elem ; 795 for j:=1:iparam do 796 sol:=sub(parm(j)=parmval(j),sol); 797 return sol 798 end ; 799 800 if sec=nonent then 801 if not fixp(omega) then 802 sol1:= for each elem in solg collect 803 begin 804 sol:=elem ; 805 for j:=1:iparam do 806 sol:=sub(parm(j)=parmval(j),sol); 807 return sol 808 end ; 809 810 if sec=entpos then 811 if fixp(omega) and (omega>0) then 812 sol1:= for each elem in solg collect 813 begin 814 sol:=elem ; 815 for j:=1:iparam do 816 sol:=sub(parm(j)=parmval(j),sol); 817 return sol 818 end ; 819 820 if sec=entneg then 821 if fixp(omega) and (omega<0) then 822 sol1:= for each elem in solg collect 823 begin 824 sol:=elem ; 825 for j:=1:iparam do 826 sol:=sub(parm(j)=parmval(j),sol); 827 return sol 828 end ; 829return sol1; 830end; 831 832procedure lectabcoef( ) ; 833%---------------------- ; 834% Lecture des coefficients de l'equation (dans l'ordre croissant des 835% derivations). 836% lecture de n : ordre de l'equation. 837% lecture des parametres (s'il apparait une variable differente de x 838% dans les coefficients). 839% les coefficients sont ranges dans la liste lcoeff (le tableau tabcoef 840% est utilise temporairement). 841% Retourne la liste { lcoeff , param } formee de la liste des 842% coefficients et de la liste des parametres (qui peut etre vide). 843% 844begin 845 scalar n, ok,iparam,lcoeff,param ; 846 847 write " " ; 848 write " ***** INTRODUCTION DES DONNEES ***** "; 849 write " " ; 850 write " L' equation est de la forme"; 851 write " a(0)(x)d^0 + a(1)(x)d^1 + .... + a(n)(x)d^n = 0 " ; 852 write " ordre de l'equation ? " ; 853 n:=xread(nil) ; 854 array tabcoef(n); 855 write " Donner les coefficients a(j)(x), j = 0..n" ; 856 for j:=0:n do tabcoef(j):=xread(nil); 857 for j:=0:n do write "a(",j,") = ",tabcoef(j); 858 write " " ; 859 write "correction ? ( oui; / non; ) " ; 860 ok:=xread(nil) ; 861 while ok eq oui do 862 << write "valeur de j :" ; j:=xread(nil) ; 863 write "expression du coefficient :";tabcoef(j):=xread(nil); 864 write "correction ?";ok:=xread(nil) ; 865 >> ; 866 867lcoeff:={tabcoef(n)}; 868for j:=n-1 step -1 until 0 do lcoeff:=tabcoef(j).lcoeff; 869if testparam(lcoeff,x) then 870<<write "nombre de parametres ? "; 871 iparam:=xread(nil); 872 if iparam neq 0 then 873 <<param:={}; 874 if iparam=1 then write "donner ce parametre :" 875 else write "donner ces parametres :"; 876 for i:=1:iparam do param:=xread(nil).param; 877 >>; 878>> else param:={}; 879clear tabcoef ; 880return {lcoeff,param}; 881end; 882 883% 884 885 886% *********************************** 887% * * 888% * UTILISATION STANDARD * 889% * * 890% *********************************** 891% 892 893 894 895procedure delire(x,k,grille,lcoeff,param) ; 896%=========================================; 897% 898% cette procedure calcule les solutions formelles d'une equation 899% differentielle lineaire homogene, a coefficients polynomiaux sur Q et 900% d'ordre quelconque, au voisinage de l'origine, point singulier 901% regulier ou irregulier ou point regulier. En fait, elle initialise 902% l'appel de la procedure NEWTON qui est une procedure recursive 903% (algorithme de NEWTON-RAMIS-MALGRANGE) 904% 905% x : variable 906% k : nombre de termes desires dans le developpement asymptotique 907% grille : les coefficients de l'operateur differentiel sont des 908% polynomes en x**grille (en general grille=1) 909% lcoeff : liste des coefficients de l'operateur differentiel (par 910% ordre croissant) 911% param : liste des parametres 912% 913% RETOURNE la liste des solutions "generales". 914% ----- 915% 916begin 917 integer prof,ordremax,ns ; 918 scalar n,l; 919 n:=length lcoeff -1; 920 array der(n),!&solution(n),!&aa(n) ; 921 array gri(20),lu(20),qx(20),equ(20),cl(20,n),clu(20,n) ; 922 array nbarete(20),xpoly(20,n),ypoly(20,n),ppoly(20,n),lpoly(20,n) ; 923 array xsq(n+1),ysq(n+1),rxm(n+1) ; 924 array ru(20,n) ,multi(20,n) ,nbracine(20) ; 925 array rac(10),ordremult(10); 926 array condprof(20),solparm(n); % liste des conditions dans Newton 927 array solequ(n); 928 on gcd ; 929 930 % initialisation du tableau cl ; 931 l:=lcoeff; 932 for i:=0:n do 933 << cl(0,i):= first l; l:=rest l;>>; 934 935 % initialisation du tableau des parametres ; 936 iparam:=length param; 937 array parm(iparam); 938 parm(0):=iparam; 939 for i:=1:iparam do parm(i):=part(param,i); 940 941 % initialisation de la grille : les coef de L sont des polynomes 942 % en x**gri(0) ; 943 gri(0):=grille ; 944 945 % substitution de d/dx par ( d/dx - (&lamb*!&u)/x**(&lamb+1) ) ; 946 der(0):=!&ff(x) ; 947 for ik:=1:n do 948 der(ik):=df(der(ik-1),x)-((!&lamb*!&u)/x**(!&lamb+1))*der(ik-1) ; 949 950 % initialisation de l'exponentielle ; 951 qx(0):=0 ; 952 953 % l'appel initial de l'algorithme NEWTON se fait avec l'operateur 954 % complet l'ordre maximum (ordremax) pour lequel on calcule le 955 % polygone NRM est n; 956 ordremax:=n ; 957 958 % initialisation de prof : prof indique le nombre d'appels recursifs 959 % de l'algorithme NEWTON ; 960 prof:=1 ; 961 962 condprof(0):={}; 963 % appel de l'algorithme NEWTON ; 964 ns:=newton(prof,ordremax,n,x,k,0) ; 965 l:=for i:=1:ns collect solequ(i); 966 clear der,!&solution,!&aa,gri,lu,qx,equ,cl,clu,nbarete,xpoly,ypoly, 967 ppoly,lpoly,xsq,ysq,rxm,tj,ru,multi,nbracine,parm ; 968 clear rac,ordremult; 969 clear condprof,solparm,solequ; 970 971 return l ; 972end; 973 974procedure testparam(l,x); 975%-----------------------; 976% l : liste des coefficients; 977% retourne t si presence de parametres (variable differente de x); 978% nil sinon; 979begin 980 scalar b,l1,l2; 981 b:=nil; l1:=l; 982 while b=nil and l1 neq{} do << l2:=coeffp({first l1},{x}); 983 for each elt in l2 do 984 <<if not numberp elt then b:=true;>>; 985 l1:=rest l1;>>; 986 return b; 987end; 988 989procedure coeffp(poly,var); 990%-------------------------; 991% poly : liste des polynomes 992% var : liste des variables 993% retourne la liste des coefficients 994begin 995 scalar l,l1 ; 996 if var={} then return poly; 997 l:={}; 998 for each elt in poly do <<l1:=coeff(elt,first var); 999 for each el1 in l1 do if el1 neq 0 then 1000 l:=append(l,{el1}) 1001 >>; 1002 return coeffp(l,rest var); 1003end; 1004 1005 1006 1007 1008procedure transformation(lcoeff,param); 1009%-------------------------------------; 1010% Entree : liste des coefficients de l'equation 1011% liste des parametres 1012% Sortie : liste des coefficients de l'equation transformee 1013begin 1014 scalar f,id,fct,fct1,coeff1,lsor; 1015 ok:=oui;coeff1:=lcoeff; 1016 while ok eq oui do <<write "derivation : 1; "; 1017 write "changement de variable : 2; "; 1018 write "changement de fonction inconnue : 3;"; 1019 write "substitution : 4;"; 1020 ichoix:=xread(nil); 1021 1022 if ichoix=1 then 1023 << write "donner le second membre : "; 1024 f:=xread(nil); 1025 write "donner le nombre de derivations : "; 1026 id:=xread(nil); 1027 coeff1:=changehom(coeff1,x,f,id) ; 1028 lsor:={coeff1,param} 1029 >>; 1030 1031 if ichoix=2 then 1032 << write "valeur de x en fonction de la nouvelle variable v ? "; 1033 fct:=xread(nil); 1034 coeff1:=changevar(coeff1,x,v,fct); 1035 coeff1:=for each elt in coeff1 collect(sub(v=x,elt)); 1036 lsor:={coeff1,param} 1037 >>; 1038 1039 if ichoix=3 then 1040 << write 1041 "valeur de y en fonction de la nouvelle fonction inconnue q ?"; 1042 fct:=xread(nil); coeff1:=changefonc(coeff1,x,q,fct); 1043 lsor:={coeff1,param} 1044 >>; 1045 1046 if ichoix=4 then 1047 << write "donner la regle de substitution , "; 1048 write "le premier membre de l'{galit{ ,puis le second : "; 1049 fct:=xread(nil); 1050 fct1:=xread(nil); 1051 lsor:=subsfonc(coeff1,param,fct,fct1); 1052 coeff1:=first lsor; 1053 >>; 1054 1055 write "transformation ? (oui;/non;) "; 1056 ok:=xread(nil); >>; 1057 return lsor; 1058end; 1059 1060procedure subsfonc(lcoeff,param,fct,fct1); 1061%----------------------------------------; 1062% Effectue la substitution de fct par fct1 1063begin 1064 scalar lsor,lsor1;integer j; 1065 lsor:= for each elt in lcoeff collect sub(fct=fct1,elt); 1066 for each elt in lsor do <<j:=j+1;write"a(",j,") = ",elt>>; 1067 lsor1:= for each elt in param do if fct neq elt then collect elt; 1068 if lsor1=0 then << 1069 write "nouvelle liste de parametres : ",{}; 1070 return {lsor,{}};>> 1071 else << 1072 write "nouvelle liste de parametres : ",lsor1; 1073 return {lsor,lsor1};>>; 1074end; 1075 1076procedure change(y,x,v,fct,exp,n); 1077%--------------------------------- 1078% exp est une expression dependant de x, de y(x), et de ses derivees 1079% d'ordre inferieur ou egal a n. 1080% change retourne la meme expression apres avoir fait le changement de 1081% variable x=fct(v). 1082begin 1083 scalar !&exp; 1084 !&hp(xt):=1/df(sub(v=xt,fct),xt); 1085 !&exp:=exp; 1086 for i:=n step -1 until 0 do !&exp:=sub(df(y,x,i)=!&d(xt,i),!&exp); 1087 !&exp:=sub(x=fct,!&exp); 1088 depend y,v; 1089 for i:=n step -1 until 0 do 1090 !&exp:=sub(df(!&fg(xt),xt,i)=df(y,v,i),!&exp); 1091 return sub(xt=v,!&exp); 1092end; 1093% 1094 1095 1096 1097% +++++++++++++++++++++++++++++++++++++++++ 1098% + + 1099% + ALGORITHME DE NEWTON + 1100% + + 1101% +++++++++++++++++++++++++++++++++++++++++ 1102%; 1103 1104 1105 1106operator !&ff,!&hp,!&fg ; 1107procedure !&d(xt,n); 1108begin 1109if n=0 then return !&fg(xt) 1110else if fixp n and (n>0) then return !&hp(xt)*df(!&d(xt,n-1),xt) ; 1111end; 1112 1113 1114procedure newton(prof,ordremax,n,x,k,ns) ; 1115%======================================= ; 1116 1117% algorithme de NEWTON-RAMIS-MALGRANGE. 1118% Cette procedure, recursive, est appelee par la procedure DELIRE. 1119% 1120% Elle NE PEUT PAS ETRE APPELEE SEULE car un certain nombre de tableaux 1121% doivent etre declares et initialises. 1122% 1123% prof : niveau de recursivite 1124% ordremax : ordre de l'operateur differentiel traite par cet appel 1125% x : variable de l'equation differentielle 1126% n : ordre de l'operateur differentiel initial 1127% k : nombre de terme du developpement asymptotique des solutions 1128% ns : nombre de solutions deja calculees lors de l'appel 1129% 1130% cette procedure retourne le nombre de solutions calculees ; 1131begin 1132 integer nba, nadep, nbsol, q ; 1133 scalar nbs,condit,sol,substitution ; 1134 1135 nbs:=ns ; 1136 1137 % construction du polygone N-R-M de l'operateur defini par 1138 % cl(prof-1,i) avec i=0..ordremax ; 1139 nba:=polygonenrm(prof,ordremax,x) ; 1140 1141 % dessin ; 1142 if lisp !*trdesir then for j:=1:nba do 1143 write xpoly(prof,j)," ",ypoly(prof,j)," ",ppoly(prof,j)," ", 1144 lpoly(prof,j) ; 1145 1146 % si la premiere arete a une pente nulle, on va calculer par FROBENIUS 1147 % lpoly(prof,1) solutions en utilisant cl(prof-1,i) ,i=0..n et 1148 % qx(prof-1) . ; 1149 % nadep = numero de la premiere arete a traiter de pente non nulle ; 1150 nadep:=1 ; 1151 1152 % si la 1ere pente est nulle : appel de frobenius et calcul des 1153 % solutions; 1154 if num(ppoly(prof,1)) = 0 then 1155 << nbsol := lpoly(prof,1) ; 1156 nouveauxaj(prof,n,x) ; 1157 condl:=condprof(prof); 1158 if lisp !*trdesir then 1159 % <<depend !&y,xt; 1160 <<write "Equation reduite : "; 1161 for i:=n step -1 until 1 do 1162 write " ",!&aa(i)," * DF(Y,XT,",i,") + "; 1163 write " ",!&aa(0)," * Y">>; 1164 % nodepend !&y,xt;>>; 1165 nbsol:=frobenius (n,xt,k) ; 1166 if nbsol neq 0 then 1167 for i:=1:nbsol do 1168 << solequ(nbs+i):={}; 1169 for each el in solparm(i) do 1170 << if length(el) > 1 then condit:=second el else condit:={}; 1171 sol:=first el; 1172 sol:=append({sub(x=xt**(1/gri(prof-1)),qx(prof-1)), 1173 gri(prof-1)},sol); 1174 solequ(nbs+i):=append (solequ(nbs+i),{{sol,condit}}); 1175 >> ; 1176 >> ; 1177 nbs:=nbs+nbsol ; 1178 nadep:=2 ; 1179 clear !&f,!°rec ; 1180 >> ; 1181 1182 % iteration sur le nombre d'aretes ; 1183 for na:=nadep:nbarete(prof) do 1184 nbs:=newtonarete(prof,na,n,x,k,nbs); 1185 % iteration sur les aretes ; 1186 1187 return nbs ; 1188end ; 1189 1190procedure newtonarete(prof,na,n,x,k,nbs); 1191%---------------------------------------; 1192begin scalar q,ordremax; 1193 q:=den(ppoly(prof,na)) ; 1194 if lisp !*trdesir then 1195 write " ",xpoly(prof,na)," ",ypoly(prof,na)," ", 1196 ppoly(prof,na)," ",lpoly(prof,na) ; 1197 1198 % calcul de la grille ; 1199 if lpoly(prof,na)=1 1200 then gri(prof):=gri(prof-1) 1201 else gri(prof):=gcd(q,1/gri(prof-1))*gri(prof-1)/q ; 1202 1203 % substitution dans l'operateur defini par cl(prof-1,i),i=0..n; 1204 lu(prof):= sub(!&lamb=ppoly(prof,na),cl(prof-1,0)*der(0)) ; 1205 for ik:=1:n do 1206 lu(prof):=lu(prof)+sub(!&lamb=ppoly(prof,na), 1207 cl(prof-1,ik)*der(ik)); 1208 1209 % decomposition de l'operateur lu ; 1210 % selon les coefficients clu(prof,i) des derivees , i=0..n ; 1211 % calcul de l'equation caracteristique ,equ(prof) : 1212 % coefficient du terme constant de clu(prof,0) ; 1213 decomplu(prof,n,x,na) ; 1214 if lisp !*trdesir then 1215 write "Equation caracteristique : ",equ(prof); 1216 1217 % calcul des racines de equ(prof) ; 1218 racinesequ(prof,na) ; 1219 1220 % iteration sur les racines de l'equation caracteristique ; 1221 for nk:=1:nbracine(prof) do 1222 << % completer l'exponentielle ; 1223 qx(prof):=qx(prof-1)+ru(prof,nk)/x**ppoly(prof,na) ; 1224 1225 % definition du nouvel operateur ; 1226 for ik:=0:n do cl(prof,ik):=sub(!&u=ru(prof,nk), 1227 clu(prof,ik)); 1228 % definition de l'ordre jusqu'auquel on calcule le nouveau 1229 % polygone-nrm : ordremax ; 1230 ordremax:=multi(prof,nk) ; 1231 1232 if lisp !*trdesir then 1233 write "Racine eq. carac. : ",ru(prof,nk); 1234 if prof <20 then nbs:=newton(prof+1,ordremax,n,x,k,nbs) 1235 else write "la profondeur 20 est atteinte :", 1236 " le calcul est arrete pour cette racine"; 1237 >> ; % fin de l'iteration sur les racines ; 1238return nbs; 1239end; 1240 1241procedure squelette (prof,ordremax,x) ; 1242%------------------------------------ ; 1243% definition du squelette du polygone de NEWTON-R-M defini par 1244% cl(prof-1,i), avec i=0..ordremax ; 1245% retourne le nombre de minima ; 1246begin 1247 scalar t00,tq,yi,cc ; 1248 integer ik,nk,nbelsq,degden,degre, rxi ; 1249 1250 condprof(prof):=condprof(prof-1); 1251 1252 % base du squelette ; 1253 % abscisse , numerotee de 1 a nbelsq ; 1254 t00:=0 ; 1255 for ik:=0 : ordremax do 1256 if cl(prof-1,ik)neq 0 then << nk:=nk+1 ; xsq(nk):=ik >> ; 1257 nbelsq:=nk ; 1258 1259 % ordonnee ; 1260 for nk:=1:nbelsq do 1261 <<tq:=sub(x=!&t**(1/gri(prof-1)),cl(prof-1,xsq(nk))) ; 1262 degden:=deg(den(tq),!&t) ; 1263 cc:=coeff(num(tq),!&t) ; 1264 ik:=0 ; 1265 while (first cc =0) do << ik:=ik+1 ; cc:= rest cc >>; 1266 ysq(nk):=(ik-degden)*gri(prof-1)-xsq(nk) ; 1267 trav1(nk):=first cc; 1268 >> ; 1269 1270 % minima successifs ; 1271 % le tableau rxm contiendra le rang de l'abscisse des minima successifs 1272 % du squelette ; 1273 % de 1 au nombre de minima ; 1274 rxm(0):=0 ; 1275 ik:=0 ; 1276 while rxm(ik)< nbelsq do 1277 <<rxi:=rxm(ik)+1 ; 1278 yi:=ysq(rxi) ; 1279 for j:=rxi+1 : nbelsq do 1280 if num(ysq(j)-yi) <= 0 then << yi:=ysq(j) ; rxi:=j >> ; 1281 ik:=ik+1 ; 1282 rxm(ik):=rxi ; 1283 >> ; 1284 return ik ; 1285end ; 1286 1287procedure polygonenrm(prof,ordremax,x) ; 1288%------------------------------------- ; 1289% construction du polygone N-R-M de l'operateur defini par cl(prof-1,i), 1290% avec i=0..ordremax ; 1291% 1292% les aretes seront numerotees de 1 a nbarete(prof) ; 1293% i=nombre d'aretes deja construites ; 1294% l'arete i est definie par : 1295% xpoly(prof,i) abscisse du sommet gauche 1296% ypoly(prof,i) ordonnee du sommet gauche 1297% ppoly(prof,i) pente de l'arete 1298% lpoly(prof,i) "longueur" de l'arete : abscisse du sommet droite - 1299% abscisse du sommet gauche; 1300% retourne le nombre d'aretes ; 1301begin 1302 scalar ydep,yfinal,pente ; 1303 integer ik,imin,jmin,nbmin,rxmin,long,xfinal,xdep,deg1,rxi ; 1304 array trav1(20); 1305 1306 nbmin:=squelette(prof,ordremax,x) ; 1307 ik:=0 ; 1308 xfinal:=xsq(rxm(1)) ; 1309 yfinal:=ysq(rxm(1)) ; 1310 xpoly(prof,1):=0 ; 1311 ypoly(prof,1):=yfinal ; 1312 ppoly(prof,1):=0 ; 1313 1314 rxi:=rxm(1); 1315 for i:=1:parm(0) do 1316 deg1:=deg1+deg(trav1(rxi),parm(i)); 1317 if deg1 neq 0 then 1318 << if lisp !*trdesir 1319 then write "Si : ",trav1(rxi)," non nul"; 1320 if (not membre({ trav1(rxi),nonnul },condprof(prof))) then 1321 condprof(prof):=cons({ trav1(rxi),nonnul },condprof(prof));>>; 1322 1323 if xfinal neq 0 then << ik:=1 ; lpoly(prof,1):=xfinal >> ; 1324 jmin:=1 ; 1325 while xfinal <ordremax do 1326 <<ik:=ik+1 ; 1327 1328 % initialisation de l'arete ik ; 1329 xpoly(prof,ik):=xfinal ; xdep:=xfinal ; 1330 ypoly(prof,ik):=yfinal ; ydep:=yfinal ; 1331 imin:=jmin+1 ; 1332 jmin:=imin ; 1333 xfinal:=xsq(rxm(imin)) ; 1334 yfinal:=ysq(rxm(imin)) ; 1335 lpoly(prof,ik):=xfinal-xdep ; 1336 ppoly(prof,ik):=(yfinal-ydep)/lpoly(prof,ik) ; 1337 1338 deg1:=0; 1339 for ii:=1:parm(0) do 1340 deg1:=deg1+deg(trav1(rxm(imin)),parm(ii)); 1341 if deg1 neq 0 then 1342 << if lisp !*trdesir 1343 then write "Si : ",trav1(rxm(imin))," non nul"; 1344 if (not membre({trav1(rxm(imin)),nonnul },condprof(prof))) then 1345 condprof(prof):=cons({ trav1(rxm(imin)),nonnul}, 1346 condprof(prof));>>; 1347 % recherche du point final de l'arete ik ; 1348 while imin < nbmin do 1349 <<imin:=imin+1 ; 1350 rxmin:=rxm(imin) ; 1351 long:=xsq(rxmin)-xdep ; 1352 pente:=(ysq(rxmin)-ydep)/long ; 1353 if num(pente-ppoly(prof,ik)) <= 0 then 1354 <<lpoly(prof,ik):=long ; 1355 ppoly(prof,ik):=pente ; 1356 xfinal:=xsq(rxmin); 1357 yfinal:=ysq(rxmin) ; 1358 jmin:=imin ; 1359 >> ; 1360 >> ; 1361 1362 >> ; 1363 nbarete(prof):=ik ; 1364clear trav1; 1365 return ik ; 1366end ; 1367 1368procedure nouveauxaj(prof,n,x) ; 1369%------------------------------ ; 1370% construction des coefficients !&aa(j) de l'operateur envoye a 1371% FROBENIUS. 1372begin 1373 scalar gr,t00,coeffs ; 1374 for i:=0:n do !&aa(i):=cl(prof-1,i) ; 1375 gr:=1/gri(prof-1); 1376 1377 % changement de x en xt**gr ; 1378 % calcul des derivees en xt ; 1379 !&hp(xt):=1/df(xt**gr,xt); 1380 1381 % calcul de l'operateur ; 1382 t00:=num( for j:=0:n sum sub(x=xt**gr,!&aa(j))*!&d(xt,j)) ; 1383 1384 % calcul des nouveaux !&aa(j) ; 1385 for j:=0:n do 1386 <<coeffs:= if j=0 then coeff(t00,!&fg(xt)) 1387 else if j=1 then coeff(t00,df(!&fg(xt),xt)) 1388 else coeff(t00,df(!&fg(xt),xt,j)); 1389 if length coeffs=1 then !&aa(j):=0 1390 else !&aa(j):=second coeffs; 1391 t00:=first coeffs 1392 >> ; 1393end ; 1394 1395procedure decomplu(prof,n,x,na) ; 1396%------------------------------- ; 1397% decomposition de l'operateur lu ; 1398% selon les coefficients clu(prof,i) des derivees , i=0..n ; 1399% calcul de l'equation caracteristique ,equ(prof) : coefficient du terme 1400% constant de clu(prof,0) ; 1401begin 1402 scalar gr,t00,tq,tj1,tj1c,coeffs ; 1403 gr:=1/gri(prof) ; 1404 t00:=num(lu(prof)) ; tq:=den(lu(prof)) ; 1405 for j:=0:n do 1406 << coeffs:= if j=0 then coeff(t00,!&ff(x)) 1407 else if j=1 then coeff(t00,df(!&ff(x),x)) 1408 else coeff(t00,df(!&ff(x),x,j)) ; 1409 if length coeffs=1 then << clu(prof,j):=0 ; t00:=first coeffs >> 1410 else << tj1:=sub(x=!&t**gr,second coeffs) ; 1411 tj1c:=coeff(tj1,!&t) ; 1412 while first tj1c =0 do tj1c:= rest tj1c; 1413 t00:=first coeffs ; 1414 if j=0 then <<clu(prof,j):=second coeffs/tq ; 1415 equ(prof):=num(first tj1c)/ 1416 !&u**(deg(num(first tj1c),!&u) 1417 -lpoly(prof,na)) 1418 >> 1419 else clu(prof,j):=second coeffs/tq ; 1420 >> ; 1421 >> ; 1422end ; 1423 1424procedure racinesequ(prof,na) ; 1425%----------------------------- ; 1426% calcul des racines de equ(prof) ; 1427begin 1428 scalar nrac ; 1429 integer nk,q1,gq,g1,dequ ; 1430 dequ:=deg(equ(prof),!&u) ; 1431 g1:=den(gri(prof-1)) ;q1:=den(ppoly(prof,na)) ; 1432 gq:=gcd(g1,q1) ; 1433 while gq > 1 do << g1:=g1/gq ;q1:=q1/gq ; 1434 gq:=gcd(g1,q1) >> ; 1435 let !&u**q1=!&t ; 1436 nrac:=racine (equ(prof),!&t) ; 1437 for ik:=1:nrac do 1438 for j:=1:q1 do 1439 << multi(prof,(ik-1)*q1+j):=ordremult(ik) ; 1440 ru(prof,(ik-1)*q1+j):=rac(ik)**(1/q1)*e**(2*(j-1)*i*pi/q1); 1441 nk:=nk+ordremult(ik) ; 1442 >> ; 1443 nbracine(prof):= nrac*q1 ; 1444 clear !&u**q1 ; 1445 if (dequ-nk) neq 0 then 1446 write "IL Y A ",ik," SOLUTIONS RELATIVES A L'ARETE " 1447 ,na," QUI NE PEUVENT PAS ETRE ATTEINTES : ", 1448 "equation a resoudre de degre >=3 " ; 1449end ; 1450 1451 1452% +++++++++++++++++++++++++++++++++++++++++ 1453% + + 1454% + ALGORITHME DE FROBENIUS + 1455% + + 1456% +++++++++++++++++++++++++++++++++++++++++ 1457%; 1458 1459 1460 1461operator !&g ; 1462% definition de !&w ; 1463% ------------------ ; 1464procedure !&w(ii,x,lambd,k); 1465if fixp k then for j:=0:k sum (df(!&g(j),lambd,ii)*x**j); 1466 1467procedure frobenius ( n,x,k ) ; 1468%============================ ; 1469 1470% Soit l'operateur differentielle : l d'ordre : n 1471% 1472% l(y)=a(n)*y(n)+a(n-1)*y(n-1)+.....a(0)*y(0) 1473% avec les a(i) = series d'ordre m en x 1474% 1475% On recherche une solution au voisinage d'un point singulier regulier 1476% de l'equation differentielle l(y)=0 sous la forme : 1477% y = x**lambda*(g(0)+g(1)*x+.....g(k)*x**k) 1478% on va determiner: 1479% - l'equation indicielle 1480% - les equations lineaires recurentes qui permettent de trouver 1481% les g(i) par identification des coefficients de x dans 1482% l'equation differentielle l(y)=0 ; 1483% 1484% Elle NE PEUT PAS ETRE APPELEE SEULE car un certain nombre de tableaux 1485% doivent etre declares et initialises. 1486% 1487% n : ordre de l'operateur 1488% x : variable 1489% k : nombre de termes du developpement asymptotique 1490% 1491% Cette procedure retourne le nombre de solutions calculees. 1492begin 1493 integer nb,nbrac,nbsolution ; 1494 scalar ss,sy, essai; 1495 equaind(n,x,k); % calcul de f(0) : equation indicielle; 1496 if lisp !*trdesir then write "Equation indicielle : ",!&f(0) ; 1497 nb:=racine (!&f(0),lambd); % calcul des racines de f(0); 1498 1499 % verification sur le calcul des racines; 1500 if nb=0 then 1501 << 1502 write "le calcul des racines est impossible dans ce cas. ", 1503 "Utilisez la version ALGEBRIQUE. "; 1504 nbsolution:=0; %cette valeur sert de test dans DELIRE; 1505 return nbsolution ; 1506 >> ; 1507 1508%etude en fonction du nombre de racines et de leur classification; 1509 nbrac:=for i:=1:nb sum ordremult(i); 1510 1511% CLASSEMENT des RACINES de l'EQUATION INDICIELLE 1512% cas particulier: 1513% ---------------- 1ou 2 racines ; 1514if nbrac=1 then 1515<< 1516 %cas d'une racine simple; 1517 nbsolution:=1; 1518 frobeniussimple(x,k,rac(1),1); 1519 solparm(1):={{{!&solution(1),rac(1)},condl} }; 1520>>; 1521 1522if nbrac=2 then << classement2r(x,k); 1523 nbsolution:=2; 1524 >> ; 1525 1526if nbrac=3 then 1527 << classement3r(x,k) ; 1528 nbsolution:=3; 1529 >>; 1530% nettoyage des variables ; 1531if nbrac>3 1532 then write "ce cas n'est pas traite. Utilisez la version ALGEBRIQUE" 1533 else for i:=0:k do clear !&g(i); 1534%fin cas ou il existe 1 ou plusieurs racines; 1535return nbsolution; 1536end ; 1537 1538procedure classement2r(x,k); 1539%--------------------------; 1540% calcul des racines lorsque l'equation indicielle a 2 racines; 1541begin 1542scalar ss,essai ; 1543if ordremult(1)=2 then rac(2):=rac(1); 1544omega:=rac(1)-rac(2); 1545if fixp(omega) then 1546 << nbsolution:=2; 1547 %if rac(2) > rac(1) then << ss:=rac(1); rac(1):=rac(2) ; 1548 % rac(2):=ss ; 1549 %modification de 10-3-93 1550 if coeffn(rac(2),i,0) > coeffn(rac(1),i,0) then << 1551 ss:=rac(1); rac(1):=rac(2); 1552 rac(2):=ss; 1553 >> ; 1554 frobeniusgeneral(x,k,nbsolution); 1555 for i:=1:nbsolution do 1556 solparm(i):={{{!&solution(i),rac(i)},condl}}; 1557 >> 1558 else 1559 if parm(0)=0 then 1560 << nbsolution:=2; 1561 frobeniussimple(x,k,rac(1),1); 1562 %pour la 2ieme solution les G(I) sont deja calcules; 1563 !&solution(2):= 1564 (for i:=0:k sum(sub(lambd=rac(2),!&g(i))*x**i)); 1565 for i:=1:nbsolution do solparm(i):={{{!&solution(i),rac(i)},condl}}; 1566 >> 1567 else 1568 << 1569%cas omega non_entier 1570 nbsolution:=2; 1571 classement2rne(x,k); 1572 >> 1573end; 1574 1575procedure classement2rne(x,k); 1576%----------------------------; 1577% calcul des racines lorsque l'equation indicielle a 2 racines et omega 1578% non-entiers 1579begin 1580scalar ss,essai ; 1581 nbsolution:=2; 1582 frobeniussimple(x,k,rac(1),1); 1583 essai:= for i:=1:k join if !&g(i)=0 then { i } else { } ; 1584 if length(essai) > 0 then essai:= ", sauf :" . essai; 1585 essai:=append({ omega, nonent }, essai); 1586 essai:=cons(essai, condl); 1587 !&solution(2):= 1588 (for i:=0:k sum(sub(lambd=rac(2),!&g(i))*x**i)); 1589 for i:=1:nbsolution do solparm(i):={{{!&solution(i),rac(i)},essai}}; 1590%cas omega >0 1591 for i:=0:k do clear !&g(i); 1592 nbsolution:=2; 1593% for i:=1:nbsolution do solparm(i):={}; 1594 frobeniusgeneral(x,k,nbsolution); 1595 essai:=cons({ omega, entpos},condl); 1596 for i:=1:nbsolution do 1597 solparm(i):=append(solparm(i),{{{!&solution(i),rac(i)},essai}}); 1598%cas omega <0 1599 for i:=0:k do clear !&g(i); 1600 nbsolution:=2; ss:=rac(1);rac(1):=rac(2);rac(2):=ss; 1601 frobeniusgeneral(x,k,nbsolution); 1602 essai:=cons({ omega, entneg},condl); 1603 for i:=1:nbsolution do 1604 solparm(i):=append(solparm(i),{{{!&solution(i),rac(i)},essai}}); 1605%cas omega =0 1606 for i:=0:k do clear !&g(i); 1607 nbsolution:=2; rac(2):=rac(1); 1608 frobeniusgeneral(x,k,nbsolution); 1609 if (not membre({ omega, entnul},condl)) then 1610 essai:=cons({ omega, entnul},condl) 1611 else essai:=condl; 1612 for i:=1:nbsolution do 1613 solparm(i):=append(solparm(i),{{{!&solution(i),rac(i)},essai}}); 1614 1615end; 1616 1617procedure classement3r(x,k) ; 1618%-------------------------- ; 1619% calcul des solutions lorsque l'equation indicielle a 3 racines ; 1620% cette procedure est appelee par FROBENIUS ; 1621begin 1622 scalar ss,sy,nbsolution ; 1623 if ordremult(1)=3 then 1624 << 1625 % cas des racines triples; 1626 rac(2):=rac(3):=rac(1) 1627 >>; 1628 if (ordremult(1)=1) and (ordremult(2)=2) 1629 then << ss:=rac(1); sy:=ordremult(1); 1630 rac(1):=rac(2); ordremult(1):=ordremult(2); 1631 rac(3):=ss; ordremult(3):=sy; 1632 >> 1633 else 1634 if ordremult(1)=2 then 1635 << 1636 %decalage de l'indice 2 puis de 1 ; 1637 rac(3):=rac(2); ordremult(3):=ordremult(2); 1638 rac(2):=rac(1); ordremult(2):=ordremult(1); 1639 >>; 1640 1641 %classement des racines ; 1642 if ordremult(1)=3 then 1643 << 1644 nbsolution:=3; 1645 frobeniusgeneral(x,k,nbsolution) 1646 >> 1647 else 1648 << % analyse des autres cas; 1649 %ordremult(1)=1; 1650 if fixp(rac(1)-rac(2)) and fixp(rac(2)-rac(3)) then 1651 << %ordonner les racines; 1652 %if rac(1)<rac(3) then << ss:=rac(1) ; 1653 % rac(1):=rac(3); rac(3):=ss ; 1654 %%modification 10-3-93 1655 !*x1:=coeffn(rac(1),i,0); 1656 !*x2:=coeffn(rac(2),i,0); 1657 !*x3:=coeffn(rac(3),i,0); 1658 if !*x1<!*x2 then << ss:=rac(1); 1659 rac(1):=rac(2); rac(2):=ss; 1660 ss:=!*x1; 1661 !*x1:=!*x2; !*x2:=ss; 1662 >> ; 1663 if !*x1<!*x3 then << ss:=rac(1); 1664 rac(1):=rac(3); rac(3):=ss; 1665 !*x3:=!*x1; 1666 >> ; 1667 if !*x2<!*x3 then << ss:=rac(2); 1668 rac(2):=rac(3); rac(3):=ss; 1669 >> ; 1670 nbsolution:=3; 1671 frobeniusgeneral(x,k,nbsolution); 1672 >>; 1673 if rac(1)=rac(2) and not fixp(rac(2)-rac(3)) then 1674 << 1675 nbsolution:=2; 1676 frobeniusgeneral(x,k,nbsolution); 1677 for i:=0:k do clear !&g(i); 1678 nbsolution:=3; 1679 frobeniussimple(x,k,rac(3),3); 1680 >>; 1681 if not fixp(rac(1)-rac(2)) and fixp(rac(2)-rac(3)) then << 1682 frobeniussimple(x,k,rac(1),3); 1683 % arranger les racines avant l'appel; 1684 rac(1):=rac(2); rac(2):=rac(3); 1685 nbsolution:=2; 1686 for i:=0:k do clear !&g(i); 1687 frobeniusgeneral(x,k,nbsolution); 1688 nbsolution:=3; 1689 >>; 1690 %cas des racines toutes distinctes n'est pas traite; 1691% if not fixp(rac(1)-rac(2)) and not fixp(rac(2)-rac(3)) then 1692 if (not fixp(rac(1)-rac(2))) and (not fixp(rac(2)-rac(3))) then 1693 << %ajout 5-5-88 ; 1694 class3rne(x,k) ; 1695 nbsolution:=3 ; 1696 >> ; 1697% write "ce cas n'est pas traite. Utilisez la version ALGEBRIQUE"; 1698 >>; 1699 for i:=1:nbsolution do 1700 solparm(i):={{{!&solution(i),rac(i)},condl}}; 1701end; 1702 1703procedure class3rne(x,k) ; 1704%----------------------- 1705begin 1706 scalar nbsolution ; 1707 if fixp(rac(1)-rac(3)) then 1708 << 1709 frobeniussimple(x,k,rac(2),3); 1710 % arranger les racines avant l'appel; 1711 rac(2):=rac(3); 1712 nbsolution:=2; 1713 for i:=0:k do clear !&g(i); 1714 frobeniusgeneral(x,k,nbsolution); 1715 nbsolution:=3; 1716 >> 1717 else 1718 << nbsolution:=3; 1719 frobeniussimple(x,k,rac(1),1); 1720 %pour la 2ieme solution les G(I) sont deja calcules; 1721 !&solution(2):= 1722 (for i:=0:k sum(sub(lambd=rac(2),!&g(i))*x**i)); 1723 !&solution(3):= 1724 (for i:=0:k sum(sub(lambd=rac(3),!&g(i))*x**i)); 1725 >>; 1726 %fin ajout; 1727end; 1728 1729procedure equaind (n,x,k) ; 1730%-------------------------- ; 1731% calcul de l'equation indicielle ; 1732% cette procedure declare un tableau f et le remplit. 1733% f(0) est l'equation indicielle ; 1734 1735% n : ordre de l'operateur 1736% x : variable 1737% k : nombre de termes demandes pour la solution ; 1738begin 1739 scalar l,denoml,ff ; 1740 integer m,di,minai,lff ; 1741 % Recherche de M=degre maximum des A(i); 1742 m:=deg(!&aa(0),x); 1743 for i:=1:n do if deg(!&aa(i),x)>m then m:=deg(!&aa(i),x); 1744 1745 array !&y(n),degrai(n),!&f(k+m+n+1); 1746 1747 % forme generale de la solution; 1748 !&y(0):=x**lambd*(for i:=0:k sum !&g(i)*x**i); 1749 1750 % determination des derivees successives de !&y; 1751 for ii:=1:n do 1752 !&y(ii):=df(!&y(ii-1),x); 1753 1754 % substitution des derivees dans l; 1755 l:=!&aa(0)*!&y(0)$ 1756 for ii:=1:n do l:=l+!&aa(ii)*!&y(ii)$ 1757 if den(l) neq 1 then << denoml:=den(l); 1758 l:=num(l); 1759 >> 1760 else denoml:=1; 1761 for ii:=0:n do 1762 << if denoml neq 1 then !&aa(ii):=!&aa(ii)*denoml; 1763 degrai(ii):= if den(!&aa(ii)) = 1 or fixp den(!&aa(ii)) 1764 then length coeff(!&aa(ii),x) -1 1765 >>; 1766 1767 % recherche du minimum entre degree(!&aa(i)) et i ; 1768 minai:=0$ 1769 maxai:=0$ 1770 1771 for ii:=0:n do 1772 << di:=degrai(ii)-ii; 1773 if (di<0) and (di<minai) then minai:=di; 1774 if (di>maxai) then maxai:=di; 1775 >>; 1776 1777 % on divise l par x**(lambd+minai); 1778 l:=l/x**(lambd+minai)$ 1779 maxai:=maxai-minai; 1780 1781 % calcul des differentes valeurs de : !&f(i); 1782 ff:=coeff(l,x)$ 1783 1784 % verification si l n'est pas divisible par : x**i; 1785 while first ff = 0 do ff:= rest ff; 1786 lff:=length ff -1; 1787 for i:=0:lff do 1788 !&f(i):=part(ff,i+1); 1789 !°rec:=maxai; 1790 1791 !&f(0):=!&f(0)/!&g(0); 1792 clear !&y,degrai ; 1793end ; 1794 1795procedure frobeniussimple (x,k,rac,nbsol) ; 1796%---------------------------------------- ; 1797% Cette procedure est particuliere a la recherche des 1798% solutions formelles d'une equation differentielle dont les solution 1799% sont simples , c.a.d. ne comportant pas de log 1800 1801% x : variable de l'equation traitee ; 1802% k : nombre de termes demande pour la solution 1803% rac : racine de l'equation indicielle 1804% nbsol : no de la solution calculee ; 1805% en fait on calcule !&solution(nbsol) ; 1806begin 1807 scalar fcoeff; array ff(k); 1808 for i:=1:k do ff(i):=!&f(i); 1809 !&g(0):=1; %choix arbitraire; 1810 for ii:=1:k do 1811 << 1812 if den ff(ii) neq 1 then ff(ii):=num(ff(ii)); 1813 if ff(ii) = 0 then !&g(ii):=0 1814 else 1815 << 1816 fcoeff:=coeff(ff(ii),!&g(ii)); 1817 !&g(ii):=-first fcoeff / second fcoeff; 1818 >>; 1819 >>; 1820 !&solution(nbsol):= (for ii:=0:k sum(sub(lambd=rac,!&g(ii))*x**ii)); 1821 clear ff; 1822end ; 1823 1824procedure frobeniusgeneral(x,k,nbsolution) ; 1825%----------------------------------------- ; 1826% x : variable de l'equation traitee ; 1827% k : nombre de termes demande pour la solution 1828% nbsolution : no de la solution calculee ; 1829begin 1830 scalar omega,fcoeff ; array ff(k); 1831 % determination des : G(i) , ce sont des fonctions de lambda ; 1832 % choix de g(0); 1833 for i:=1:k do ff(i):=!&f(i); 1834 if nbsolution = 2 then 1835 << 1836 if rac(1)=rac(2) then !&g(0):=1 1837 else 1838 << 1839 % on suppose que les racines sont ordonnees de facon croissante 1840 % c.a.d. rac(1)>rac(2); 1841 omega:=rac(1)-rac(2); 1842 !&g(0):=sub(lambd=lambd+omega,!&f(0)); 1843 >>; 1844>>; 1845if nbsolution = 3 then 1846<< 1847 omega:=rac(1)-rac(3); 1848%? if omega<0 then omega :=-omega; 1849% probleme pour la determination de G(0) - A revoir et verifier; 1850 1851 !&g(0):=for i:=1:omega product( sub(lambd=lambd+i,!&f(0)) ); 1852>>; 1853 1854for i:=1:k do 1855<< 1856 %rappel K fixe (nombre de terme demande); 1857 ff(i):=num(ff(i)); 1858 if ff(i) = 0 then !&g(i):=0 1859 else 1860 << 1861 fcoeff:=coeff(ff(i),!&g(i)); 1862 !&g(i):=-first(fcoeff)/second(fcoeff); 1863 >>; 1864>>; 1865 1866if lisp !*trdesir then 1867 << write "Solution en l'indeterminee lambda : "; 1868 factor x; 1869 write !&w(0,x,lambd,k); 1870 remfac x>>; 1871 1872%determination des solutions; 1873 1874if rac(1)=rac(2) then 1875<< 1876 !&solution(1):=sub(lambd=rac(1),!&w(0,x,lambd,k)); 1877 !&solution(2):=sub(lambd=rac(1),!&w(0,x,lambd,k)) 1878 *log(x) 1879 + sub(lambd=rac(1),!&w(1,x,lambd,k)); 1880>> 1881else 1882 << 1883 !&solution(1):=sub(lambd=rac(1),!&w(0,x,lambd,k)); 1884 if parm(0)=0 then 1885 !&solution(2):=sub(lambd=rac(2),!&w(0,x,lambd,k)) 1886 *log(x) + 1887 sub(lambd=rac(2),!&w(1,x,lambd,k)) 1888 else 1889 !&solution(2):=!&w(0,x,lambd,k) 1890 *log(x) + !&w(1,x,lambd,k); 1891 >>; 1892 1893 if nbsolution = 3 then 1894 !&solution(3):=sub(lambd=rac(3),!&w(0,x,lambd,k)) 1895 *(log x)**2 1896 + 2*sub(lambd=rac(3),!&w(1,x,lambd,k)) 1897 *log(x) 1898 + sub(lambd=rac(3),!&w(2,x,lambd,k) ) ; 1899 clear ff; 1900end ; 1901 1902 1903 1904% +++++++++++++++++++++++++++++++++++++++++++++++ 1905% + + 1906% + PROCEDURES UTILITAIRES + 1907% + + 1908% +++++++++++++++++++++++++++++++++++++++++++++++ 1909%; 1910 1911 1912 1913procedure racine(f,x) ; 1914%-------------------- ; 1915% procedure qui calcule les racines quelconques ( et leur ordre de 1916% multiplicite ) d'une equation algebrique ; 1917% 1918% f : on cherche les racines de l'equation algebrique f(x)=0 1919% x : variable 1920% 1921% rac : tableau a une dimension des racines distinctes (de 1 a nbrac) 1922% ordremult : tableau correspondand de leur ordre de multiplicite 1923% cette procedure retourne le nombre de racines distinctes ; 1924begin 1925 integer nbrac ; 1926 scalar sol, multsol ; 1927 nbrac:=0 ; 1928 sol:=solve(f,x); 1929 multsol:=multiplicities!* ; 1930 for each elt in sol do 1931 if lhs(elt) = x then 1932 << nbrac:=nbrac+1 ; 1933 ordremult(nbrac):=first(multsol); 1934 multsol:=rest(multsol) ; 1935 rac(nbrac):=rhs(elt) ; 1936 >> 1937 else multsol:=rest(multsol) ; 1938 return nbrac ; 1939end ; 1940 1941procedure membre(elt,list); 1942%------------------------; 1943begin 1944scalar bool; 1945for each w in list do if w=elt then bool:= t; 1946return bool; 1947end; 1948 1949symbolic ; 1950if !*trdesir then << 1951terpri() ; terpri() ; 1952princ " DESIR : solutions formelles d'equations differentielles" ; 1953terpri() ; 1954princ " lineaires homogenes au voisinage de zero, point " ; 1955terpri() ; 1956princ " singulier regulier ou irregulier, ou point regulier" ; 1957terpri() ; terpri() ; 1958princ " Version 3.3 - Novembre 1993 " ; terpri() ; 1959princ " Appel par desir(); " ; terpri() ; 1960terpri() ; 1961>>; 1962 1963algebraic ; 1964on gcd ; 1965 1966endmodule; 1967 1968end; 1969