(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)(* DDMF is distributed under CeCILL-B license. *)(* The service LocalSolutionsOrigin computes all transseries solutions *)(* of a differential equation at 0. A transseries solution consists of *)(* a Puiseux series of log-power type, multiplied by some exponential *)(* factor and a power factor x^alpha with alpha an algebraic number. *)(* Here is an overview how several services call each other to compute *)(* the local solutions of a differential equation at some point: *)(* LocalSolutionsPoint: transseries solutions at an arbitrary point *)(* -> LocalSolutionsOrigin: transseries solutions at the origin *)(* -> LogSeriesAlpha: log-power series with x^alpha factor *)(* -> LogSeriesZero: log-power series solutions *)INCLUDE"preamble.ml"lettitle (_, _, _, _, _) = <:text<LocalSolutionsat theOrigin>>(* Compute the generalized exponents. This proc basically calls *)(* DEtools[gen_exp] and postprocesses its output a little bit, *)(* by simplifying RootOf expressions, extracting the ramification *)(* indices, and flattening the lists. The final output is of the *)(* form [[e1, r1], [e2, r2], ...] where ej is a generalized *)(* exponent (i.e., exp(int(ej/x,x)) is an exponential factor) and *)(* rj is the corresponding ramification index (an integer). *)letgeneralized_exponents eqn1 y1 x1 = << (proc(eqn, y, x, $) local genexp, i, g, a, r, s,T, pf, ram; genexp :=DEtools[gen_exp](eqn, y(x),T); # get ridofRootOf-s that have no index parameter genexp := map(allvalues, genexp);fori from 1tonops(genexp)dog := genexp[i]; # we expect the last elementtobeofthe form a*T^r=x.ifnot(match(op(1, g[-1]) = a*T^r,T,'s')andevalb(op(0, g[-1]) = `=`)andevalb(op(2, g[-1]) = x))thenerror("Generalized exponent different than expected.")endif; # ram is the ramification indexandpf the above prefactor a. ram := subs(s, r); pf := subs(s, a); g := subsop(-1 =NULL, g); g := map(a -> subs(T= x/a, g), [allvalues(RootOf(_Z^ram - pf))]); g := map(a -> [a, ram],ListTools:-Flatten(g)); genexp[i] := g;enddo; returnListTools:-FlattenOnce(genexp);endproc)($(eqn1), $(y1), $(x1)) >>(* This function is called from the main service in order to loop over *)(* all exponential factors. It returns a pair (doc, sols) that consists *)(* of the documented trace of the calculations and an Ocaml list *)(* containing the solutions corresponding to this exponential factor. *)letdeal_with_exp_factor parallel eqn y x u n exp_factor ram =(* Perform the transformation that corresponds to the ramification index. *)lettrans_eqn = << collect( evala(gfun:-algebraicsubs($(eqn), $(y) - $(x)^$(ram), $(y)($(x)))), {diff, $(y)}, factor) >>inletram_text =if<:bool< $(ram) = 1 >>then<:par<>>else<:par<The$(t_ent:Glossary.g"ramification index")forthis exponential factor is <:isymb<$(ram)>>.Therefore, removing the exponential factorasit is would resultina differential equationwithhigher orderandproduce spurious solutions.Inordertoavoid this the substitution <:imath< $(symb: x) \mapsto $(symb: x)^$(symb: ram)>> is performed, which leadstothe differential equation <:dmath< $(symb: trans_eqn) = 0. >> >>in(* Remove the exponential factor from the differential equation. *)lettrans_eqn = << collect( evala(gfun:-`diffeq*diffeq`( $(trans_eqn), gfun:-holexprtodiffeq(1/$(exp_factor), $(y)($(x))), $(y)($(x)))), {diff, $(y)}, factor) >>inletexp_text =if<:bool< $(exp_factor) = 1 >>then<:par<>>else<:par<Theexponential factor <:isymb< $(exp_factor) >> is removed by applying the transformation <:imath< << $(y)($(x)) >> \mapsto << simplify(1/$(exp_factor)) * $(y)($(x)) >> >>.Thisresultsinthe differential equation <:dmath< $(symb: trans_eqn) = 0. >> >>in(* Compute all solutions of the form x^alpha * logseries. *)letseries_text, lsers =ifparallelthenDC.inline_service (LogSeriesAlpha.descr (trans_eqn, y, x, u, n)None),LogSeriesAlpha.obj ((trans_eqn, y, x, u, n), ())elseLogSeriesAlpha.doc_obj trans_eqn y x u ninletlsers =CommonTools.symb_list_of_symb lsersin(* Include the exponential factor into the log-power-series solutions, *)(* and undo the ramification substitution. *)List.iter (funa -> <:unit< $(a):-exp_factor := subs($(x) = $(x)^(1/$(ram)), $(exp_factor)); $(a):-alpha := $(a):-alpha / $(ram); $(a):-finite_part := subs($(x) = $(x)^(1/$(ram)), $(a):-finite_part); $(a):-ramification := $(ram); >>) lsers; (ram_text @:@ exp_text) @@@ series_text, lsers(* Compute all transseries solutions of a differential equation at 0. *)(* The functions u[i](n) are used to encode the coefficient sequences *)(* in the series expansions. *)letdoc_obj eqn y x u n =letgenexp = generalized_exponents eqn y xin(* Determine equivalence classes (same exponential) among the *)(* generalized exponents. *)letgenexp = <<DDMF:-equivalence_classes( $(genexp), (a, b) -> not(has(simplify(a[1] - b[1]), $(x)))andevalb(a[2] = b[2])) >>in(* Now take some representative for each equivalence class. *)(* The sort is only to enforce deterministic behaviour of maple. *)letgenexp = << sort(map(a -> op(1, a), $(genexp))) >>in(* Delete the constant part of each polynomial; this corresponds to a *)(* power of x only and should go into the alpha of LogSeriesAlpha. *)(* Of course, one could use this information directly, but we decided to *)(* perform the extraction of alpha ourselves using the indicial equation. *)letgenexp = << map(a -> [a[1] - subs($(x) = infinity, a[1]), a[2]], $(genexp)) >>in(* Transform the Laurent polynomials into exponentials. *)letexp_factors =CommonTools.symb_list_of_symb << map(a -> exp(int(a[2] * a[1] / $(x), $(x))), $(genexp)) >>in(* Extract the ramification indices. *)letrams =CommonTools.symb_list_of_symb << map(a -> a[2], $(genexp)) >>in(* For each exponential factor compute the log-powerseries solutions. *)letparallel =List.length exp_factors > 1inletfcall = deal_with_exp_factor parallel eqn y x u ninletcontents, sol_list =List.split (List.map2 fcall exp_factors rams)in(* The exponential factors as they appear in the result. *)letexp_factors =List.map (funa -> << $(List.hd a):-exp_factor >>) sol_listinletnexps =List.length exp_factorsin(* If there are nontrivial exponential factors, write something about *)(* how they were obtained. *)letintro =ifnexps = 1 && <:bool< $(List.hd exp_factors) = 1 >>thenDC.ent_nullelse<:par<Themodular algorithm describedin$(t_ent:Bibliography.b"ClHo04") is employedtodetermine that the local solutions involve the exponential factor$(str:Wording.ending_of_int nexps) >> @:@ (Wording.enumeration_of_maples exp_factors) @:@ <:par<.>>in(* If there are more than just one exponential factor, create subsections. *)letcontent =ifnexps = 1thenList.hd contentselseList.fold_left (@@@)DC.ent_null (List.map2 (funa b ->DC.section <:text<Solutionswithexponential factor <:isymb< $(a) >>>> b) exp_factors contents)inintro @@@ content,CommonTools.symb_of_symb_list (List.flatten sol_list) let_serviceLocalSolutionsOrigin(eqn : diffeq maple) (y : name maple) (x : name maple) (u : name maple) (n : name maple) :DC.sec_entities * any maplewith{ title = title } = doc_obj eqn y x u n

Generated by GNU Enscript 1.6.5.90.