# LogSeriesZero.ml

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

(* The service LogSeriesZero computes all log-power series solutions   *)
(* of a differential equation at 0, i.e., solutions of the form        *)
(*   \sum_{n=0}^\infty \sum_{j=0}^{t-1} u_j(n) \ln(x)^j x^n.           *)
(* For this purpose it is assumed that at least one such solution      *)
(* exists, i.e., that the origin is not an irregular singular point,   *)
(* and that the equation has been transformed such that the smallest   *)
(* integer root of the indicial equation is 0.                         *)

(* 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<Log-Power Series at the Origin>>

(* This function is called from the main service in order to loop over  *)
(* all canonical solutions. Each canonical solution is characterized by *)
(* case = << [k, d] >> where k is the index of the root                 *)
(* n_0 = int_roots[k] and d is an index labeling the different          *)
(* solutions that correspond to this root (i.e., whose series expansion *)
(* starts at n_0).                                                      *)
(* It returns a pair (doc, sol) that consists of the documented trace   *)
(* of the calculations and the solution as a Maple Record.              *)
let deal_with_case x u n bigQ q int_roots mults case =
let k = << $(case)[1] >> and d = <<$(case)[2] >> in
let n0 = << $(int_roots)[$(k)] >>
and m = << $(mults)[$(k)] >> in

(* Compute t early for a mathematically sound local expansion. *)
(* (The origin of t is explained later.) *)
let t =
<< $(d) + 1 + add($(mults)[i], i = $(k) + 1 .. nops($(int_roots))) >> in
let solution_form =
<:par<
<:imath< << ln($(x))^$(d) * $(x)^$(n0) >> +
O(<< ln($(x))^($(t) - 1) * $(x)^($(n0) + 1) >>) >>\
>> in

let intro =
<:par<Since the indicial equation has the integer root <:isymb< $(n0) >> >> @:@ (if <:bool<$(m) = 1 >>
then <:par<, there exists a $(t_ent:Glossary.g "formal logarithmic sum") >> else <:par< with multiplicity <:isymb<$(m) >>, there exist <:isymb< $(m) >> linearly independent$(t_ent:Glossary.g
~display:(Some "formal logarithmic sums")
"formal logarithmic sum")
>>)
@:@
<:par<
whose lowest exponent of <:isymb< $(x) >> is exactly <:isymb<$(n0) >>.
>> @:@
(if <:bool< $(m) = 1 >> then <:par<>> else <:par<This subsection deals with the >> @:@ (Wording.ordinal <:int<$(d) + 1 >>) @:@
<:par< of these solutions, which takes the form >> @:@
solution_form @:@ <:par<.>>) in

(* Compute t and describe what's going on. *)
let m1 =
CommonTools.symb_list_of_symb
<< [seq($(mults)[i], i =$(k)+1 .. nops($(mults)))] >> in let sum_mults = List.map (fun a -> "+" ^ (Maple.serialization_of_t a)) m1 in let sum = List.fold_left (^) (Maple.serialization_of_t d) sum_mults in let t_text = if <:bool<$(t) = 1 >> then <:par<>> else
<:par<
There, the highest possible exponent of <:isymb< ln($(x)) >> is determined by the sum of the current index (which is <:isymb<$(d) >>)
and the multiplicities of all integer roots of the indicial
equation that are greater than the current root, <:isymb< $(n0) >>; this gives <:imath<$(str: sum) = << $(t) - 1 >> >>. >> in (* This is basically formula (3.9), *) (* where u[k](n) has been substituted by k!*u[k](n-n0) *) let s = <<$(n) = $(n) +$(n0) >> in
(* Since the Q's contain the parameters of the given special function, *)
(* use of local variables is essential here! *)
let recsys = <<
proc(bigQ, u, t, s, n) local i, j, l;
[seq(
factor(subs(s, bigQ[1,1])) * u[i](n) =
add(factor(subs(s, -bigQ[1,j-i+1]) * binomial(j, i)) *
u[j](n), j = i+1 .. t-1) +
u[j](n-l), j = i .. t-1), l = 1 .. $(q)), i = 0 .. t-1)] end proc($(bigQ), $(u),$(t), $(s),$(n)) >> in

(* Explain and display the recurrences. *)
let recs1 = CommonTools.symb_list_of_symb recsys in
let recs2 =
List.map
(fun a -> <:par<<:dmath< $(symb: a), >>>>) (List.rev (List.tl recs1)) in let recs3 = recs2 @ [ <:par<<:dmath< <<$(List.hd recs1) >>.>>>> ] in
let recs = List.fold_left (@:@) <:par<>> recs3 in
let smnd = << add($(u)[i]($(n))*ln($(x))^i, i = 0 ..$(t)-1) * $(x)^$(n) >> in
let ansatz =
if <:bool< $(n0) = 0 >> then <:par<<:dmath< \sum_{$(symb: n)=0}^\infty $(symb: smnd) >>>> else <:par<<:dmath< <<$(x)^$(n0) >> \sum_{$(symb: n)=0}^\infty $(symb: smnd) >>>> in let recsys_text = <:par<Applying the differential equation to the \$(t_ent:Glossary.g "ansatz")>>
@:@ ansatz @:@
<:par<and comparing equal powers of <:isymb< $(x) >> >> @:@ (if <:bool< nops($(recsys)) = 1 >>
then <:par<leads to the recurrence >>
else <:par<and <:isymb< ln($(x)) >> leads to the system of recurrences >>) @:@ recs @:@ <:par<An alternative expression for >> @:@ (if <:bool< nops($(recsys)) = 1 >>
then <:par<this recurrence>>
else <:par<these recurrences>>) @:@
<:par< \
is in terms of the polynomials <:imath< Q_j($(symb: n)) >> computed above and reads <:dmath< \sum_{\ell=0}^$(symb: q) \sum_{j=i}^<< $(t)-1 >> {j\choose i} Q_{\ell}^{(j-i)}(<<$(n) + $(n0) >>) \,$(symb: u)_j($(symb: n)-\ell) = 0 \quad (0 \leq i <$(symb: t)).>>
>> in

(* Initial values *)
let inits = <<
[
seq($(u)[j](0) = 0, j = 0 ..$(d)-1),
$(u)[$(d)](0) = 1,
seq($(u)[j](0) = 0, j =$(d)+1 .. $(t)-1) ] >> in let inits_text1 = <:par<Since the solution should be of the form >> @:@ solution_form @:@ <:par<\ , the initial value$(str: Wording.ending_of_int <:int< $(t) >>) for <:imath<$(symb: n) = 0 >> must be chosen as follows:
>> @:@
(Wording.enumeration_of_maplelist inits) @:@
<:par<. >> in
let arb_values =
<< [seq(seq($(u)[j]($(int_roots)[i] - $(n0)), j = 0 ..$(mults)[i]-1),
i = $(k)+1 .. nops($(int_roots))) ] >> in
let narb = <:int< nops($(arb_values)) >> in let inits_text2 = if narb = 0 then <:par<>> else <:par<A closer inspection of the recurrences above reveals that >> @:@ (Wording.enumeration_of_maplelist arb_values) @:@ <:par< never occur>> @:@ (if narb = 1 then <:par<s. Its value>> else <:par<. Their values>>) @:@ <:par< therefore may be chosen arbitrarily, <:imath<0>> for example. >> in let inits = << [op($(inits)), op(map(a -> a = 0, $(arb_values)))] >> in (* Compute the remaining initial values using the recurrences and *) (* the assumption that the u[i]'s at negative arguments are zero. *) let eqns = << [seq(op(subs($(n) = i, $(recsys))), i = 1 .. max($(q)-1, max($(int_roots)) -$(n0)))] >> in
let eqns = << subs($(inits),$(eqns)) >> in
let eqns = << subs(
[seq(seq($(u)[j](i) = 0, j = 0 ..$(t)-1), i = 1-$(q) .. -1)],$(eqns)) >> in
let vars =
<< indets($(eqns), specfunc(integer, {seq($(u)[j], j = 0 .. $(t)-1)})) >> in let sols = << solve($(eqns), [op($(vars))])[1] >> in let is = << {seq(i, i = 0 ..$(q)-1), op(map(a -> a - $(n0),$(int_roots)))} >> in
let sols = << select(a -> member(op(1, op(1, a)), $(is)),$(sols)) >> in
let inits = << [op($(inits)), op($(sols))] >> in

(* If we give additional initial values, display them and say where they *)
(* come from. *)
let u0orj =
if <:bool< $(t) = 1 >> then <<$(u)[0]($(n)) >> else <<$(u)[j]($(n)) >> in let inits_text3 = if <:bool< nops($(sols)) = 0 >> then <:par<>> else
<:par<
But for the sake of completeness and convenience, all initial values
<:isymb< $(symb: u0orj) >>, >> @:@ (if <:bool<$(t) = 1 >>
then <:par<>>
else <:par<<:imath< 0 \leq j \leq << $(t) - 1 >>>> and >>) @:@ <:par<<:imath< 0 \leq <<$(n) >> < $(symb: q) >> must be given>> @:@ (if <:bool<$(t) = 1 >> then <:par<. >> else
<:par<,
as well as those <:isymb< $(symb: u0orj) >> for which the leading coefficient of the respective recurrence vanishes. >>) @:@ <:par< The missing values, however, can be obtained from the previous ones using the recurrence$(str: Wording.ending_of_seq recsys)
and the assumption that
<:imath< $(symb: u0orj) = 0 >> for <:imath<$(symb: n) < 0 >>.
This produces <:dmath< << op($(sols)) >>. >> >> in (* Discard sequences whose initial values are all 0. *) let oldt = t and inits, recsys, t = ref inits, ref recsys, ref t in while <:bool< {0, op(map(a -> op(2, a), select(has,$(!inits), $(u)[$(!t)-1])))} = {0}
>>
do
inits := << remove(has, $(!inits),$(u)[$(!t)-1]) >>; recsys := << eval(subsop(-1 = NULL,$(!recsys)), $(u)[$(!t)-1] = 0) >>;
t := << $(!t) - 1 >> done; let inits, recsys, t = !inits, !recsys, !t in let zero_seq_text = if <:bool<$(t) = $(oldt) >> then <:par<>> else <:par<It turns out that all initial values of the sequence>> @:@ <:par<$(str: Wording.ending_of_int <:int< $(oldt) -$(t) >>) >> @:@
(Wording.enumeration_of_maplelist
<< [seq($(u)[i]($(n)), i = $(t) ..$(oldt)-1)] >>) @:@
<:par< are zero. Hence >> @:@
(if <:bool< $(oldt) -$(t) = 1 >>
then <:par<this sequence is>>
else <:par<these sequences are>>)
@:@
<:par< identically zero and can be discarded. >> in

let solution = <<
Record(
'exp_factor' = 1,
'alpha' = $(n0), 'finite_part' = 0, 'closed_form' = undefined, 'variable' =$(x),
'ramification' = 1,
'log_powers' = [seq(i, i = 0 .. $(t)-1)], 'sequence_name' =$(u),
'index_name' = $(n), 'recurrences' = ListTools:-Reverse($(recsys)),
'initial_conditions' = sort($(inits))) >> in intro @:@ t_text @:@ recsys_text @:@ inits_text1 @:@ inits_text2 @:@ inits_text3 @:@ zero_seq_text, solution (* Compute all power series solutions of eqn at the origin of the form *) (* \sum_{n=0}^\infty \sum_{j=0}^{t-1} u_j(n) \ln(x)^j x^n *) (* where it is assumed that 0 is the smallest integer root *) (* of the indicial equation. *) (* The names u and n have to be given as arguments. *) (* Equation numbers and notation are according to *) (* Joris van der Hoeven: Fast Evaluation of Holonomic Functions Near *) (* and in Regular Singularities (JSC 31, 717-743, 2001). *) let doc_obj eqn y x u n = let ord = << DDMF:-ode_order($(eqn)) >> in

(* Transform eqn into Euler operator notation delta = x*D_x. *)
(* In view of the correspondence (delta(f))_n = n*u(n) we use delta = n. *)
let eqn_euler =
<< subs([$(y)($(x)) = 1, seq(diff($(y)($(x)), $(x)$ j) =
$(x)^(-j) * add(combinat[stirling1](j, i) *$(n)^i, i = 1..j),
j = 1 .. $(ord))],$(eqn)) >> in
let eqn_euler =
<< expand($(x)^(-ldegree(expand($(eqn_euler)), $(x))) *$(eqn_euler)) >> in
let eu = << theta[$(x)] >> in let euler_text = <:par< It is obtained by rewriting the differential equation into operator notation in terms of the Euler operator <:imath< <<$(eu)(f($(x))) >> =$(symb: x)f'($(symb: x)) >>: <:dmath< << collect(subs($(n) = $(eu),$(eqn_euler)), $(eu), factor) >>.>> >> in (* Compute the polynomials Q that are introduced in Formula (3.7) *) (* and their derivatives: Q[i+1,j+1] corresponds to Q_i^{(j)}. *) (* Q_0 = Q[1,1] is the indicial polynomial. *) let bigQ = << [seq(collect(subs($(n) = $(n) - i, coeff($(eqn_euler), $(x), i)),$(n)),
i = 0..degree($(eqn_euler),$(x)))] >> in
let bigQ =
<< map(a -> [seq(diff(a, [$(n)$ j]), j = 0 .. degree($(bigQ)[1],$(n)))],
$(bigQ)) >> in let q = << nops($(bigQ)) - 1 >> in
let qs =
List.map
(fun a ->
<:par<<:dmath<
Q_$(int: a-1)($(symb: n)) = << $(bigQ)[$(int: a), 1] >>,
>>>>)
(CommonTools.range <:int< $(q) + 1 >>) in let bigQ_text = <:par< From this form, it can be easily read off how the differential equation acts on the Taylor coefficients of <:imath< f($(symb: x)) = \sum_{$(symb: n)=0}^\infty f_$(symb: n) $(symb: x)^$(symb: n) >>.
The rules
<:dmath< [$(symb:x)^$(symb:n)] \, ($(symb:eu) f($(symb:x))) =
$(symb:n) f_$(symb:n), \quad
[$(symb:x)^$(symb:n)] \, ($(symb:x) f($(symb:x))) = f_{$(symb:n)-1} >> imply that the <:isymb<$(n) >>th Taylor coefficient of the series
obtained by substituting <:imath< f($(symb: x)) >> into the differential equation is given by >> @:@ (match <:int<$(q) >> with
| 0 -> <:par<<:dmath< Q_0($(symb: n)) f_$(symb: n) >>>>
| 1 -> <:par<<:dmath< Q_0($(symb: n)) f_$(symb: n) +
Q_1($(symb: n)) f_{$(symb: n)-1} >>>>
| 2 -> <:par<<:dmath< Q_0($(symb: n)) f_$(symb: n) +
Q_1($(symb: n)) f_{$(symb: n)-1} +
Q_2($(symb: n)) f_{$(symb: n)-2} >>>>
| _ -> <:par<<:dmath< Q_0($(symb: n)) f_$(symb: n) + \dots +
Q_$(symb: q)($(symb: n)) f_{$(symb: n)-$(symb: q)} >>>>)
@:@
<:par<
for some polynomials <:imath< Q_j($(symb: n)) >>, <:imath< 0 \leq j \leq$(symb: q) >>.
(For the notation <:imath< [$(symb:x)^$(symb:n)] \, (\dots) >>, see
$(t_ent:Glossary.g "formal power series").) Each <:imath< Q_j($(symb: n)) >>
is obtained as the coefficient of <:imath< $(symb: x)^j >> in the differential operator above, where <:isymb<$(eu) >> has been replaced
by <:imath< $(symb: n) - j >>, leading to: >> @:@ (List.fold_left (@:@) <:par<>> qs) @:@ <:par< and it is seen in this way \ that the indicial polynomial is <:imath< Q_0 >>. >> in (* The integer roots of the indicial equation and their multiplicities. *) let all_roots = << [solve($(bigQ)[1,1] = 0, $(n))] >> in let int_roots = << sort(select(type,$(all_roots), integer)) >> in
let int_roots = << DDMF:-equivalence_classes($(int_roots), =) >> in let mults = << map(nops,$(int_roots)) >>
and int_roots = << map(a -> op(1, a), $(int_roots)) >> and sorted_roots = << convert(sort(convert($(all_roots), set)), list) >> in

(* Display the indicial equation and its root(s). *)
let ind_text =
<:par<
The indicial equation is
<:dmath< << $(bigQ)[1,1] >> = 0 , >> whose >> @:@ (if <:bool< nops($(sorted_roots)) = 1 >>
then <:par<only root is >>
else <:par<roots are >>)
@:@ (Wording.enumeration_of_maplelist sorted_roots) @:@ <:par<.>> in

(* All pairs [k, d] of indices for which a canonical solution f^{[n0,d]} *)
(* exists where n0 = int_roots[k]. *)
let all_cases = CommonTools.symb_list_of_symb
<<
[seq(seq([k, d], d = 0 .. $(mults)[k] - 1), k = 1 .. nops($(int_roots)))]
>>
in

let contents, sols =
List.split
(List.map (deal_with_case x u n bigQ q int_roots mults) all_cases) in

(* Create subsections only if there are several solutions. *)
let sec_titles =
List.map
(fun a -> <:text<The canonical solution
<:imath< f^{[ << $(int_roots)[$(a)[1]] >>, << \$(a)[2] >> ]} >>>>)
all_cases in
let content =
if List.length sols = 1
then List.hd contents
else List.fold_left
(@@@)
DC.ent_null
(List.map2 DC.section sec_titles contents) in

(ind_text @:@ euler_text @:@ bigQ_text) @@@ content,
CommonTools.symb_of_symb_list sols

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