# LocalSolutionsOrigin.ml

```(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)

(* 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"

let title (_, _, _, _, _) = <:text<Local Solutions at the Origin>>

(* 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).       *)
let generalized_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 rid of RootOf-s that have no index parameter
genexp := map(allvalues, genexp);

for i from 1 to nops(genexp) do
g := genexp[i];

# we expect the last element to be of the form a*T^r=x.
if not(match(op(1, g[-1]) = a*T^r, T, 's')
and evalb(op(0, g[-1]) = `=`)
and evalb(op(2, g[-1]) = x))
then error("Generalized exponent different than expected.")
end if;

# ram is the ramification index and pf 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;
end do;
return ListTools:-FlattenOnce(genexp);
end proc)(\$(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.   *)
let deal_with_exp_factor parallel eqn y x u n exp_factor ram =

(* Perform the transformation that corresponds to the ramification index. *)
let trans_eqn =
<< collect(
evala(gfun:-algebraicsubs(\$(eqn), \$(y) - \$(x)^\$(ram), \$(y)(\$(x)))),
{diff, \$(y)}, factor) >> in
let ram_text = if <:bool< \$(ram) = 1 >> then <:par<>> else
<:par<
The \$(t_ent: Glossary.g "ramification index") for this exponential
factor is <:isymb<\$(ram)>>. Therefore, removing the exponential
factor as it is would result in a differential equation with higher
order and produce spurious solutions. In order to avoid this
the substitution
<:imath< \$(symb: x) \mapsto \$(symb: x)^\$(symb: ram)>>
is performed, which leads to the differential equation
<:dmath< \$(symb: trans_eqn) = 0. >>
>> in

(* Remove the exponential factor from the differential equation. *)
let trans_eqn =
<< collect(
evala(gfun:-`diffeq*diffeq`(
\$(trans_eqn),
gfun:-holexprtodiffeq(1/\$(exp_factor), \$(y)(\$(x))), \$(y)(\$(x)))),
{diff, \$(y)}, factor) >> in
let exp_text = if <:bool< \$(exp_factor) = 1 >> then <:par<>> else
<:par<
The exponential factor <:isymb< \$(exp_factor) >> is removed
by applying the transformation
<:imath< << \$(y)(\$(x)) >> \mapsto
<< simplify(1/\$(exp_factor)) * \$(y)(\$(x)) >> >>.
This results in the differential equation
<:dmath< \$(symb: trans_eqn) = 0. >>
>> in

(* Compute all solutions of the form x^alpha * logseries. *)
let series_text, lsers =
if parallel then
DC.inline_service (LogSeriesAlpha.descr (trans_eqn, y, x, u, n) None),
LogSeriesAlpha.obj ((trans_eqn, y, x, u, n), ())
else
LogSeriesAlpha.doc_obj trans_eqn y x u n
in
let lsers = CommonTools.symb_list_of_symb lsers in

(* Include the exponential factor into the log-power-series solutions, *)
(* and undo the ramification substitution. *)
List.iter
(fun a ->
<: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. *)
let doc_obj eqn y x u n =

let genexp = generalized_exponents eqn y x in

(* Determine equivalence classes (same exponential) among the *)
(* generalized exponents. *)
let genexp =
<<
DDMF:-equivalence_classes(
\$(genexp),
(a, b) -> not(has(simplify(a[1] - b[1]), \$(x))) and evalb(a[2] = b[2]))
>> in

(* Now take some representative for each equivalence class. *)
(* The sort is only to enforce deterministic behaviour of maple. *)
let genexp = << 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. *)
let genexp =
<< map(a -> [a[1] - subs(\$(x) = infinity, a[1]), a[2]], \$(genexp)) >> in

(* Transform the Laurent polynomials into exponentials. *)
let exp_factors =
CommonTools.symb_list_of_symb
<< map(a -> exp(int(a[2] * a[1] / \$(x), \$(x))), \$(genexp)) >> in

(* Extract the ramification indices. *)
let rams = CommonTools.symb_list_of_symb << map(a -> a[2], \$(genexp)) >> in

(* For each exponential factor compute the log-powerseries solutions. *)
let parallel = List.length exp_factors > 1 in
let fcall = deal_with_exp_factor parallel eqn y x u n in
let contents, sol_list = List.split (List.map2 fcall exp_factors rams) in

(* The exponential factors as they appear in the result. *)
let exp_factors =
List.map (fun a -> << \$(List.hd a):-exp_factor >>) sol_list in
let nexps = List.length exp_factors in

(* If there are nontrivial exponential factors, write something about *)
(* how they were obtained. *)
let intro =
if nexps = 1 && <:bool< \$(List.hd exp_factors) = 1 >>
then DC.ent_null
else
<:par<
The modular algorithm described in \$(t_ent:Bibliography.b "ClHo04")
is employed to determine 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. *)
let content =
if nexps = 1
then List.hd contents
else List.fold_left (@@@) DC.ent_null
(List.map2
(fun a b ->
DC.section
<:text<Solutions with exponential factor <:isymb< \$(a) >>>> b)
exp_factors
contents) in

intro @@@ content, CommonTools.symb_of_symb_list (List.flatten sol_list)

let_service LocalSolutionsOrigin
(eqn : diffeq maple)
(y : name maple)
(x : name maple)
(u : name maple)
(n : name maple) :
DC.sec_entities * any maple with { title = title } =
doc_obj eqn y x u n
```

Generated by GNU Enscript 1.6.5.90.