# LogSeriesAlpha.ml

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

(* The service LogSeriesAlpha computes all log-power series solutions  *)
(* of a differential equation at 0 including possibly an x^alpha       *)
(* factor, i.e., solutions of the form        *)
(*   x^alpha * \sum_{n=0}^\infty \sum_{j=0}^{t-1} u_j(n) \ln(x)^j x^n, *)
(* where alpha is an algebraic number.                                 *)
(* It is assumed that at least one such solution exists, i.e., that 0  *)
(* is not an irregular singular point of the differential equation.    *)

(* 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<Series Solutions with Power Factor at the Origin>>

(* This function is called from the main service in order to loop over  *)
(* all shifts that are determined by the indicial equation. It returns  *)
(* a pair (doc, sols) that consists of the documented trace of the      *)
(* calculations and an Ocaml list sols containing the solutions that    *)
(* correspond to the given shift.                                       *)
let deal_with_shift parallel eqn y x u n shift =

(* Remove the power factor. *)
let trans_eqn =
<< collect(
evala(gfun:-diffeq*diffeq(
$(eqn),$(x) * diff($(y)($(x)), $(x)) +$(shift) * $(y)($(x)), $(y)($(x)))),
{diff, $(y)}, factor) >> in let fac_text = if <:bool<$(shift) = 0 >> then <:par<>> else
<:par<
The power factor <:isymb< $(x)^$(shift) >> is removed by means of the
transformation
<:imath< << $(y)($(x)) >> \mapsto << $(x)^$(shift) * $(y)($(x)) >> >>
which results in the differential equation
<:dmath< $(symb: trans_eqn) = 0. >> >> in (* Compute the log-power-series solutions. *) let content, psers = if parallel then DC.inline_service (LogSeriesZero.descr (trans_eqn, y, x, u, n) None), LogSeriesZero.obj ((trans_eqn, y, x, u, n), ()) else LogSeriesZero.doc_obj trans_eqn y x u n in let psers = CommonTools.symb_list_of_symb psers in (* Include the power factor into the log-power-series solutions. *) List.iter (fun a -> <:unit<$(a):-alpha := $(shift) +$(a):-alpha >>) psers;

fac_text @@@ content, psers

(* Compute all local solutions of a differential equation at 0, which *)
(* are of the type                                                    *)
(*   x^\alpha \sum_{n=0}^\infty \sum_{j=0}^{t-1} u_j(n) \ln(x)^j x^n  *)
(* where alpha is an algebraic number.                                *)
(* The names u and n have to be given as arguments.                   *)
let doc_obj eqn y x u n =

let intro =
<:par<
Studying the $(t_ent: Glossary.g "indicial equation") of the differential equation above explains what exponents of <:isymb<$(x) >> can be involved in local series solutions.
>> in

(* Compute the roots of the indicial equation and group them into *)
(* classes determined by integer shift equivalence.               *)
(* The equivalence classes are sorted by their first elements,    *)
(* particularly for enforcing deterministic behaviour of maple.   *)
let ind_eqn = << evala(gfun:-indicialeq(
gfun:-formatdiffeq([$(eqn),$(y)($(x))]),$(x), $(n))) >> in let ind_roots = << [solve($(ind_eqn) = 0, $(n))] >> in let ecs = << DDMF:-equivalence_classes($(ind_roots), (a, b) -> type(a - b, integer)) >> in
let ecs = << map(sort, $(ecs)) >> in let ecs = CommonTools.symb_list_of_symb << sort($(ecs), (a, b) -> [a[1], b[1]] = sort([a[1], b[1]])) >> in
let necs = List.length ecs in

(* Create the corresponding content. *)
let ind_eqn_text =
if necs = 1 && <:bool< $(List.hd ecs)[1] = 0 >> then (* The indicial equation is always displayed by the LogSeriesZero. *) (* In this case the same equation would be displayed twice in a row. *) <:par<>> else <:par< The indicial equation in this case is <:dmath<$(symb: ind_eqn) = 0 >> and its
>> @:@
(if <:bool< nops($(ind_roots)) = 1 >> then <:par<only root is >> else <:par<roots are >>) @:@ (Wording.enumeration_of_maples (CommonTools.symb_list_of_symb ind_roots)) @:@ <:par<.>> in let ecs_text = if necs <= 1 then <:par<>> else <:par< These roots are grouped into classes determined by integer shift equivalence, i.e., all elements in an equivalence class differ by an integer only. Thus the equivalence classes are >> @:@ (Wording.enumeration_of_maples ecs) @:@ <:par<, each of which is treated separately in one of the following subsections. >> in (* Extract the "smallest" element s from each equivalence class, so that *) (* all others are of the form s+n with n being a positive integer. *) let shifts = List.map (fun a -> <<$(a)[1] + min(map(b -> b - $(a)[1],$(a))) >>) ecs in

let parallel = List.length shifts > 1 in
let fcall = deal_with_shift parallel eqn y x u n in
let contents, sols = List.split (List.map fcall shifts) in

(* Create subsections if not all indicial roots are are *)
(* in the same equivalence class. *)
let content =
if necs = 1
then List.hd contents
else List.fold_left (@@@) DC.ent_null
(List.map2
(fun a b ->
DC.section
<:text<Solutions with power factor <:isymb< $(x)^$(a) >>>> b)
shifts
contents) in

(intro @:@ ind_eqn_text @:@ ecs_text) @@@ content,
CommonTools.symb_of_symb_list (List.flatten sols)

let_service LogSeriesAlpha
(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.