DiffeqToRecurrence.ml

```(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)
(* however this file is distributed under the LGPL version 2.1. *)

INCLUDE "preamble.ml"

let title (_, _, _, _, notation) =
<:text<Finding a Recurrence for the Power Series Coefficients
of <:imath< \$(str: notation) >>>>

(* The following exception is used from goodinitvalues_rec to diffeqtorec. *)
exception No_solution

(* The following exception is used internally in goodinitvalues_rec only. *)
(* FC: This is how I understood the former code when porting it to ocaml. *)
(* But it is really strange. *)
exception Return_without_display of any maple

(*
# `goodinitvalues/rec`
# Input:  a recurrence in the format returned by gfun:-formatrec
#    the unknown sequence and its variable
#    some initial conditions
#    a boolean flag cleanup
#    (optional) an integer p
#    (optional) maxsing = nonnegint
# Output: a set of equalities u(k)=v_k, from which all the other values
#    (smaller than maxsing) can be deduced by solving the recurrence for its
#     maximal index.
#    This set is continued up to the pth term if p is given.
#    New variables _C[i] are introduced when a value is arbitrary.
#    When the flag is true, the unnecessary entries u(k)=_C[i] are removed
#    The result is an ERROR when no initial condition can be found.
*)

let goodinitvalues_rec recurrence u n ini flag argsix =

(* :=proc (recurrence, u, n, ini, flag, arg6)
option `Copyright (c) 1992-2009 by Algorithms Project, INRIA France. \
local n0, order, i, inds, minind, maxind, sys, r, sol, b, a, j, k,
rej, maxsingularity, gb, termorder, dorej, inds2; *)

<:unit<
# TODO: remove the following, when these are local vars
unassign('n0', 'i', 'inds', 'minind', 'maxind', 'sys', 'r', 'sol', 'b',
'a', 'j', 'k',
'rej', 'maxsingularity', 'gb', 'termorder', 'dorej', 'inds2');

recOrder:=nops(\$(recurrence))-2;  # SG: replaced order by recOrder
maxind:=recOrder-1;
if type(\$(ini),'set') then
inds:=map(op,indets(\$(ini),\$(u)('integer')))
else
inds:={}; maxind:=max(maxind,\$(ini))
end if;
maxsingularity:=NULL;

if \$(argsix)<>NULL then    # SG: replaced "nargs=6"
if type(\$(argsix),integer) then
maxind:=max(maxind,\$(argsix))
elif type(\$(argsix),identical('maxsing')=nonnegint) then
maxsingularity:=op(2,\$(argsix))
else error "invalid argument",\$(argsix)
fi
end if;

# The initial conditions are u(minind)..\$(u)(maxind).

if maxsingularity<>NULL then
n0 := gfun:-firstnonzero(
subs(\$(n) = \$(n)-recOrder,
\$(recurrence)[nops(\$(recurrence))]),
\$(n),
maxsingularity);
else
n0 := gfun:-firstnonzero(
subs(\$(n) = \$(n)-recOrder,
\$(recurrence)[nops(\$(recurrence))]),
\$(n));
end if;

# u(n0-1) cannot be deduced from the previous ones.
maxind:=max(maxind,op(inds),n0-1);
minind:=min(op(inds),0);
r:=gfun:-makerec(\$(recurrence),\$(u),\$(n));
>> ;

(* SG: breakFlag instead of "break" statement in the following loop *)
<:unit< breakFlag := false: >> ;

(* Some cases will return a computed value but no do any display. *)
(* An exception Return_without_display is used in this case. *)
(* The 'catch' is at the very bottom of the function. *)
try

let sol_found = ref false and j = ref 1 in
while not !sol_found && !j < 3 do
(*
<* j := \$(int:j): *> (* added by SG; j is used below *)
*)

if !j = 1 then
<:unit<
sys := { op(\$(ini)),
seq(subs(\$(n)=i,r), i=minind..maxind-recOrder) } :
>>
else
<:unit<
sys := { op(\$(ini)),
seq(subs(\$(n)=i-recOrder,r),
i={\$minind+recOrder..maxind} minus inds) } :
>> ;

if <:bool< sys = {} >> then
raise
(Return_without_display
(if not flag then
<< { seq(\$(u)(i) = _C[i], i=minind..maxind) } >>
else
<<
{ seq(\$(u)(i) = _C[i],
i={\$minind..maxind} minus {\$0..recOrder-1}) }
>>)) ;

<:unit<
a := gfun:-systomatrix(sys,[seq(\$(u)(i),i=minind..maxind)],'b');
>> ;
sol_found :=
(* FIXME: We need a DynaMoW:-Raise badly enough! *)
(* But will extraction of Maple code be easy? *)
<:bool<
try
sol := LinearAlgebra:-LinearSolve(a,Vector['column'](b))
catch "inconsistent system" :
sol := NULL
end try ;
evalb(sol <> NULL)
>> ;

if not !sol_found && !j = 2 then
<:unit< error "no valid initial conditions": >> ;

incr j
done ;

(* SG: added the following if *)
if not !sol_found then
raise No_solution ;

(* added by SG; TODO should become local var *)
<:unit< sysOld :=sys : >> ;

(* TODO: replaced k by kk in the following; should become local var *)
<:unit<
sol:=convert(sol,list);
#    inds:=indets(sol,_t[anything]);
#    inds:=indets(sol,'typeindex'( 'anything', 'suffixed'('_t', 'integer') ));
#  SG: added the Or; before, e.g., _t[1] would not be replaced
inds := indets(sol,
'typeindex'('anything',
Or(identical('_t'),
'suffixed'('_t', 'integer')) ));

# replace the _t[anything] by _C[anything] depending on flag
dorej:=\$(bool:flag) and evalb(j=1); # SG: added evalb

j:=max(op(map(op,indets([\$(recurrence),\$(ini)],
_C['anything']))));

if j=-infinity then j:=-1 end if;
for i in inds do
if member(i,sol,'kk') and (not dorej or nops(select(has,sol,i))>1
or kk<minind-1 or kk+minind>recOrder) then
j:=j+1;
sol:=subs(i=_C[minind+j],sol);
rej[i]:=NULL
else rej[i]:=kk
end if
end do;

sys:={seq(\$(u)(i+minind-1)=sol[i],i={\$1..nops(sol)} minus
{seq(rej[i],i=inds)})};
>> ;

(* The end of this function is structured as: *)
(*   begin looong_display end, begin looong_symbolic end *)
(* catch Return_without_display symbolic_system -> *)
(*  <:par<>>, symbolic_system *)
(* The unusual begin/end are used to mitigate the overload of nested *)
(* ifs and structures. *)

(* Start of long text computation, to be returned. *)
begin
if <:bool< n0 > recOrder >> then
<:par<The first index from which on the leading coefficient is non-zero
is <:isymb< \$(n) = n0-recOrder >>.
The value <:isymb< \$(u)(n0-1) >>
cannot be deduced from the previous ones by the recurrence. Hence >> @:@
(let bound = <:int< n0-recOrder-1 >> in
if <:bool< n0-recOrder-1 = 0 >> then
<:par<the instance <:imath< \$(symb:n) = 0 >> of the recurrence becomes
an initial condition.>>
else
<:par<the instances
<:imath< \$(symb:n) = 0, \dots,~\$(symb:n) = \$(int:bound) >>
of the recurrence become inital conditions.>>) @:@
(if <:bool< maxind > minind >> then
<:par<The initial sequence elements
<:imath< \$(symb:u)(minind), \dots,~\$(symb:u)(maxind) >>
thus satisfy the >> @:@
(if <:bool< nops(sysOld) > 1 >> then
<:par<system of linear equations>>
else
<:par<linear equation>>) @:@
<:par<<:dsymb<
op({ op(map(x->x=0, select(x-> not type(x,'equation'), sysOld))),
op(select(x->type(x,'equation'), sysOld)) })
>>>>
else
(if <:bool< maxind = minind >> then
<:par<The sequence element <:imath< \$(symb:u)(minind) >> thus
satisfies
<:dsymb<
op({ op(map(x->x=0, select(x-> not type(x,'equation'), sysOld))),
op(select(x->type(x,'equation'), sysOld)) })
>>
>>
else
<:par<>>))
else
(if <:bool< sys <> sysOld >> then
<:par<Because of the constraints on the initial conditions, >> @:@
if <:bool< maxind > minind >> then
<:par<the sequence elements
<:imath< \$(symb:u)(<< minind >>), \dots,~\$(symb:u)(<< maxind >>) >>
satisfies the >> @:@
(if <:bool< nops(sysOld) > 1 >> then
<:par<system of linear equations>>
else
<:par<linear equation>>) @:@
<:par<<:dmath<
<< op({ op(map(x->x=0, select(x-> not type(x,'equation'), sysOld))),
op(select(x->type(x,'equation'), sysOld)) }) >> .
>>>>
else
(if <:bool< maxind = minind >> then
<:par<the sequence element <:imath< \$(symb:u)(minind) >> satisfies
<:dsymb<
op({ op(map(x->x=0,
select(x-> not type(x,'equation'), sysOld))),
op(select(x->type(x,'equation'), sysOld)) })
>>
>>
else
<:par<>>)
else
<:par<>>) @:@

(if <:bool< sys <> sysOld >> then
(if <:bool< nops(sysOld) > 1 >> then
<:par<The solution of this system is >>
else
<:par<The solution of this equation is >>) @:@
<:par<<:dmath< << op(sys) >> . >>>>
else
<:par<>>)
end,

(* Clean the _C[] of initial conditions of the type _C[1] +/- _C[2]. *)
begin
if <:bool<
hastype(remove(type,sys,\$(u)('anything')='name'),_C['anything'])
# limitations due to the cost of Groebner basis computation:
and max(op(map(degree,[seq(op(2,i),i=sys)],
indets(sys,_C['anything']))))<3 and not has(sys,'RootOf')
# also it has to be a system of polynomials
and type([seq(op(2,i),i=sys)],
'list'('polynom'('rational',indets(sys,_C['anything']))))
>> then
let _ =
<:unit<
inds:=indets(sys,_C['anything']);
inds2:=sort([op(map(op,indets(sys,\$(u)('anything'))))]);
sys:=subs([seq(\$(u)(i)=\$(u)[i],i=inds2)],
{seq(op(1,i)-op(2,i),i=sys)}):
# find algebraic relations between the u[i]
termorder:=lexdeg([op(inds)],[seq(\$(u)[i],i=inds2)]);
gb:=remove(hastype,Groebner:-Basis(sys,termorder),_C['anything']);
sol:=subs(solve({op(gb)},{seq(\$(u)[i],
i=inds2[Groebner:-HilbertDimension(
gb,termorder)-nops(inds)..-1])}),
[seq(\$(u)[i],i=inds2)]);
j:=-1;
for i to nops(sol) do
if sol[i]=\$(u)[inds2[i]] then
if dorej and nops(select(has,sol,sol[i]))=1 then
rej[i]:=i
else
rej[i]:=NULL;
j:=j+1;
sol:=subs(\$(u)[inds2[i]]=_C[j],sol);
end if
end if
end do;
>> in
<<
{seq(\$(u)(inds2[i])=sol[i],i={\$1..nops(sol)} minus
{seq(rej[i],i=1..nops(sol))})}
>>
else
<< sys >>
end

(* Deals with the cases that require/permit no display. *)
with Return_without_display symbolic_system ->
<:par<>>, symbolic_system

let diffeqtorec_doit big_r y z u k iniconds yDisplay zDisplay =
(* := proc(big_r,y,z,u,k,iniconds)
option `Copyright (c) 1992-2009 by Algorithms Project, INRIA France. \
local l, ini, i, rec, j, minordrec, maxordrec, m, r, dr1, inhdeg, inhpart,
p, rr, cont; *)

(* SG: added parameters yDisplay, zDisplay. They are used in DynaMoW
instead of y and z, as diffeqtorec always sets the latter to Y and Z. *)

(* TODO local variables *)
(*  annotate Maple <* local nn, firstk ; *> ; *)

<:unit<
if has(\$(big_r),\$(k)) then
error \$(k),`cannot appear in the differential equation`
elif has(\$(big_r),\$(u)) then
error \$(u),`cannot appear in the differential equation` end if;

# initial conditions
# iniconds is there. SG: originally the test was "nargs=6"
if \$(iniconds)<>NULL then
l:={seq(op(2,op(0,op(0,i))),
i=indets([\$(iniconds),\$(big_r)],
`gfun/initeq`(\$(y)))
minus {\$(y)(0),'D'(\$(y))(0)})};
ini:=[\$(y)(0)=\$(u)(0),
'D(\$(y))(0)'=\$(u)(1),
seq(`@@`(D,i)(\$(y))(0)=\$(u)(i)*i!,i=l)];
r:=subs(ini,\$(big_r)); ini:=subs(ini,\$(iniconds))
else r:=\$(big_r); ini:={} end if;
>> ;

let par1 =
if <:bool< ini <> {} >> then
(if <:bool< nops(iniconds) = 1 >> then
<:par<The initial condition of the differential equation implies>>
else
<:par<The initial conditions of the differential equation imply>>) @:@
<:par<<:dmath< << op(ini) >> . >>>>
else
<:par<>>
in

<:unit<
# In very special cases, this loop makes it possible to return
# an inhomogeneous equation of lower order.
# Ex: z*(-1+z)^3*(D@@2)(y)(z)+(-1+z)^3*D(y)(z)-(-1+z)^3*y(z)-z*(z-3)
inhdeg := 0;
if r[1]<>0 then
for inhdeg from 0 do
p:=1-\$(z);
for i from 2 to nops(r) while degree(p,\$(z))=1 do
p:=gcd(p,r[i])
end do;
if not has(p,\$(z)) then break end if;
r:=[r[1],seq(quo(r[i],1-\$(z),\$(z)),i=2..nops(r))]
end do
end if;
# main loop
## this fixes a problem with the change in the definition of degree(0)
## in V.5
rr:=map(proc(x) if x=0 then 1 else x end if end,r);
>> ;

<:unit<
eqnWithPowerSeries :=
DDMF:-PowerSeries(\$(k), i, mul(k-j,j=0..i-1)*\$(u)(k),
\$(zDisplay), -i),
i=0..nops(r)-2) =
subs(\$(z)=\$(zDisplay),
-r[1])/(1-\$(zDisplay))^inhdeg
>> ;

let inhom_deg = <:int< inhdeg >> in

let par2 =
(if inhom_deg > 0 then
<:par<Dividing by the factor
<:imath< (1-\$(symb:zDisplay))^\$(int:inhom_deg) >> and replacing
>>
else
<:par<Replacing >>) @:@
<:par<<:imath< \$(symb:yDisplay)(\$(symb:zDisplay)) >> by the power series
<:isymb< DDMF:-PowerSeries(\$(k),0,\$(u)(\$(k)),\$(zDisplay),0) >>
in the differential equation results in
<:dmath< << eqnWithPowerSeries >> . >>
>>
in

<:unit<
minordrec:=min(seq(i-degree(op(i+2,rr),\$(z)),i=0..nops(r)-2));
maxordrec:=max(seq(i-ldegree(op(i+2,rr),\$(z)),i=0..nops(r)-2));
rec:=Array(minordrec..maxordrec,'storage'='sparse');
for i from 2 to nops(r) do
for j from ldegree(op(i,r),\$(z)) to degree(op(i,r),\$(z)) do
rec[i-2-j]:=rec[i-2-j]+coeff(op(i,r),\$(z),j)*
mul(k+m,m=1-j..i-2-j) end do end do;

# SG: "expand" removed from main loop, put below
# SG: added remPoly (polynomial built from initial power series terms,
#      which do net reflect the full recurrence)

maxDeg := max(seq(degree(op(i,r)), i=2..nops(r)));
kInit := max(0,max(seq(degree(op(i,r))-(i-2), i=2..nops(r))));
mul(kk+i-j-nn, nn=0..i-1)*\$(z)^kk, kk=j..kInit-1), j=0..maxDeg-1),
i=0..nops(r)-2),\$(z));

# inhomogeneous part of the differential equation
if r[1]=0 then dr1:=-1; inhpart:=0
else
dr1:=degree(r[1],\$(z));
if inhdeg<>0 then # r[1] stands for r[1]/(1-z)^inhdeg
j=1-i..inhdeg-i-1),
i=0..dr1))/(inhdeg-1)!;
r := subsop(1=series(r[1]/(1-\$(z))^inhdeg,\$(z),
max(dr1,maxordrec-minordrec)+1),r)
else
inhpart:=0
end if
end if;
>> ;

(* We now collect powers and shift summation indices. An equation is
displayed only if it differs from the preceding one. *)

let par3 =
if inhom_deg > 0 then
(let _ =
(* TODO intermediate step on the rhs *)
<:unit<
eqnCollectPowers :=
subs(\$(k)=\$(k)-i,rec[i])*\$(u)(k),
\$(zDisplay), -i),
i=minordrec..maxordrec) +
subs(\$(z)=\$(zDisplay), remPoly) =
- DDMF:-PowerSeries(\$(k), 0, inhpart, \$(zDisplay), 0) :
>> in
if <:bool< eqnCollectPowers <> eqnWithPowerSeries >> then
<:par<We expand the right-hand side, using the binomial theorem.
After collecting powers of <:isymb< zDisplay >>
on the left-hand side, the equation becomes
<:dmath< << eqnCollectPowers >> . >>
>>
else
<:par<>>) @:@
(let _ =
<:unit<
eqnShift :=
rec[i]*\$(u)(k+i), \$(zDisplay), 0),
i=minordrec..maxordrec) +
subs(\$(z)=\$(zDisplay), remPoly) =
-DDMF:-PowerSeries(\$(k), 0, inhpart, \$(zDisplay), 0):
>> in
if <:bool< eqnShift <> eqnCollectPowers >> then
<:par<Now we shift summation indices, so that the same
power <:isymb< \$(zDisplay)^\$(k) >>
appears in all sums, and we obtain
<:dsymb< eqnShift >>
>>
else
<:par<>>)
else
(let _ =
<:unit<
eqnCollectPowers :=
subs(\$(k)=\$(k)-i, rec[i]) * \$(u)(k),
\$(zDisplay), -i), i=minordrec..maxordrec) +
subs(\$(z)=\$(zDisplay), remPoly) =
subs(\$(z)=\$(zDisplay), -r[1]):
>> in
if <:bool< eqnCollectPowers <> eqnWithPowerSeries >> then
<:par<After collecting powers of <:imath< \$(symb:zDisplay) >>,
the equation becomes
<:dmath< << eqnCollectPowers >> . >>
>>
else
<:par<>>) @:@
(let _ =
<:unit<
eqnShift :=
rec[i]*\$(u)(k+i), \$(zDisplay), 0),
i=minordrec..maxordrec) +
subs(\$(z)=\$(zDisplay), remPoly) =
subs(\$(z)=\$(zDisplay), -r[1]):
>> in
if <:bool< eqnShift <> eqnCollectPowers >> then
<:par<Now we shift summation indices, so that the same
power <:isymb< \$(zDisplay)^\$(k) >>
appears in all sums, and we obtain
<:dmath< << eqnShift >> . >>
>>
else
<:par<>>)
in

<:unit<
for i from minordrec to maxordrec do rec[i] := expand(rec[i]) end do;

# SG: Changed "nargs=6" to "iniconds<>NULL".
if ini<>{} or \$(iniconds)<>NULL or dr1<>-1 then
seq(rec[j]*\$(u)(i+j),j=max(minordrec,-i)..maxordrec)]),
i=0..dr1)],`+`))};
for i from dr1+1 while i<-minordrec or subs(k=i,{seq(rec[maxordrec-j],
j=0..min(i-dr1-1,maxordrec-minordrec))})={0} do
j=max(minordrec,-i)..maxordrec),inhpart]),`+`)}
end do;
end if;
>> ;

let par4 =
let sOrNo = Wording.ending_of_seq << addini >> in
<:par<Comparing coefficients of <:isymb< \$(zDisplay)^\$(k) >>
gives the >> @:@
(if <:bool< addini <> {} >> then
<:par<initial value constraint\$(str:sOrNo)
<:dsymb< op(map(x -> x=0, addini)) >>
and the >>
else
<:par<>>) @:@
(let _ =
(* TODO: recurrence holds for smaller values, if inhdeg>0 *)
(* firstk added by SG; the first k, for which the recurrence holds *)
<:unit< firstk := max(kInit, dr1 + 1) >> in
<:par<following recurrence, valid for <:isymb< \$(k) >= firstk >>:
<:dmath<
<< add(rec[i]*\$(u)(\$(k)+i), i=minordrec..maxordrec) + inhpart = 0 >> .
>>
>>)
in

<:unit<
rec := gfun:-listprimpart(
subs(k=k-minordrec, [inhpart,seq(rec[i],i=minordrec..maxordrec)]),
\$(k),
'cont');
>> ;

let par5 =
if <:bool<
evalb(minordrec <> 0 or (cont <> 1 and not type(cont,'symbol')))
>>
then
(if <:bool<
evalb(minordrec <> 0 and cont <> 1 and not type(cont,'symbol'))
>>
then
<:par<Cancelling the common
factor <:isymb< subs(\$(k)=\$(k)+minordrec, cont) >>
and shifting <:isymb<\$(k)>> by <:isymb< -minordrec >>,>>
else
if <:bool< cont <> 1 and not type(cont,'symbol') >> then
<:par<Cancelling the common
factor <:isymb< subs(\$(k)=\$(k)+minordrec,cont) >>,>>
else
<:par<Shifting <:isymb<\$(k)>> by <:isymb< -minordrec >>,>>) @:@
(let _ =
<:unit< firstk := max(kInit, dr1+1) + minordrec: >> in
<:par< we find that, for <:isymb< \$(k) >= firstk >>,
<:dmath<
<< add(rec[nn]*\$(u)(\$(k)+nn-2), nn=2..maxordrec-minordrec+2) = 0 >> .
>>
>>)
else
<:par<>>
in

<:unit<
# SG: changed "nargs=6" to "iniconds<>NULL"
if has(cont,\$(k)) and (ini<>{} or \$(iniconds)<>NULL) then
l:=-1;
for j in select(type,gfun:-myisolve(cont,\$(k)),nonnegint) do
if j>=i then # this i is the last one in the previous loop
l:=l+1;
fi
od
fi;
>> ;

let par6 =
if <:bool< addini <> {} >> then
<:par<Due to the zeros of the canceled factor,
the following sequence values are arbitrary:
<:dmath< << addini >> . >>
>>
else
<:par<>>
in

<:unit<
while nops(rec)>2 and rec[nops(rec)]=0 do
rec:=subsop(nops(rec)=NULL,rec) end do;
>> ;

(par1 @@@ par2 @@@ par3 @@@ par4 @@@ par5 @@@ par6,
(* firstk added by SG *)
<< map(collect, rec, \$(k)), ini, cont, firstk >>)

(* diffeqtorec
Input:  eqn: differential equation (for example output of algeqtodiffeq)
y(z): its unknown function
u(n): the name of the sequence of Taylor coefficients at the origin
ini: (optional) boolean indicating whether initial conditions should
be computed
content_rec: (optional) will be assigned the content that was removed from
the recurrence. This is useful when the diff. eqn. is only the homogenous
part of the actual one.
Output: the linear recurrence satisfied by u(n) *)
(* Added by SG: the second member of the output ist the first n for which *)
(* the recurrence holds. *)

let diffeqtorec eqn yofz uofk ini =

(* TODO: {ini::boolean:=false,return_content::boolean:=false,
content_rec::name:=dummy} *)
(*  diffeqtorec:=proc (eqn,yofz,uofk, {ini::boolean:=false,
return_content::boolean:=false,content_rec::name:=dummy}) *)
(*   option `Copyright (c) 1992-2009 by Algorithms Project, INRIA France. \
(*   local iniconds, f, y, z, u, k, Y, Z, rec, inirec, cont; *)

(* retVal added by SG *)
(*  symbolic retVal ; *)

(* TODO remove, as soon as iniconds is local *)
<:unit< unassign('iniconds'): >> ;

let inirec = ref << "empty" >> in

<:unit<
# This also checks that >=3 args are passed.
gfun:-getname(\$(uofk),'u','k'):
>> ;

if <:bool< \$(bool:ini) or type(\$(eqn),'set') >> then
<:unit< f := gfun:-formatdiffeq([\$(eqn),\$(yofz)],'y','z','iniconds'): >>
else
<:unit< f := gfun:-formatdiffeq([\$(eqn),\$(yofz)],'y','z'): >> ;

let par1 =

<:par<The function under consideration satisfies the differential equation
<:dsymb< add(op(i+2,f)*diff(y(z),[z\$i]), i=0..nops(f)-2) = -op(1,f) >>
>> @:@

(let sOrNo = Wording.ending_of_seq << iniconds >> in
if <:bool< type(iniconds, 'set') and not iniconds = {} >> then
<:par<with initial condition\$(str:sOrNo)
<:dmath< << op(iniconds) >> . >>
>>
else
<:par<>>) @:@
<:par<Our goal is to find a recurrence for the Taylor
coefficients <:isymb< \$(uofk) >>
of <:isymb< \$(yofz) >>.>>

in

(* Avoid problems when u=y or u=z, ... *)
<:unit< f := subs([y=Y,z=Z], f): >> ;

if <:bool< \$(bool:ini) or type(\$(eqn),'set') >> then
<:unit< iniconds := subs(y=Y, iniconds): >>
else
<:unit< iniconds := NULL: >> ;

(* SG: added d2recReturn; TODO should become local *)
let par2, d2recReturn =
diffeqtorec_doit
<< f >>
<< Y >>
<< Z >>
<< u >>
<< k >>
<< iniconds >>
<< y >>
<< z >> in

<:unit< rec := op(1,[\$(d2recReturn)]): >> ;
inirec := << op(2,[\$(d2recReturn)]) >> ;
<:unit<
cont := op(3,[\$(d2recReturn)]):
firstk := op(4,[\$(d2recReturn)]):
>> ;

let par3 =
if <:bool< iniconds<>{} or iniconds<>NULL: >> then
let p, ir =
goodinitvalues_rec
<< rec >>
<< u >>
<< k >>
!inirec
true
<< NULL >> in
begin
inirec := ir ;
p
end
else
<:par<>>
in

(* SG: Reassign u and k; this is a work-around. So far there is no way *)
(* to have local variables, hence k may have been erroneously modified *)
(* by a procedure call. *)
<:unit< gfun:-getname(\$(uofk),'u','k'): >> ;

let result, par4 =
(* if added by SG *)
(* FIXME: FC, I changed the coding-style so that not a false is returned *)
(* but an exception is raised. *)
let the_inirec = !inirec in
if <:bool< \$(the_inirec) = false >> then
<< false, infinity >>,
<:par<The initial conditions are contradictory, hence
the differential equation above has no power series >> @:@
(if <:bool< type(iniconds,'set') and iniconds<>{} >> then
<:par<solution with the given initial values.>>
else
<:par<solution.>>)
else
let r = << gfun:-makerec(rec,u,k,\$(the_inirec)), firstk >> in
r,
if <:bool< type(op(1, [\$(r)]), 'set') and nops(op(1,[\$(r)])) > 1 >>
then
let the_inirec = << firstk >> in
<:par<Summarizing, the coefficients <:imath< u(k) >>
of any power series solution of the differential equation above
satisfy, for <:imath< k \geq \$(symb:the_inirec) >>,
the recurrence relation
<:dsymb< op(1,op(1,[\$(r)])) = 0 >>
>> @:@
(* TODO: Give references to the steps where the initial conditions *)
(* occur. *)
(let sOrNo =
Wording.ending_of_int (<:int< nops(op(1,[\$(r)])) >> - 1) in
<:par<with initial condition\$(str:sOrNo)
<:dmath< << op(2..nops(op(1,[\$(r)])), op(1,[\$(r)])) >> . >>
>>)
else <:par<>> in
(result, par1 @@@ par2 @@@ par3 @@@ par4)

let_service DiffeqToRecurrence
(eqn : diffeq maple)
(yofz : fcall maple)
(uofk : fcall maple)
(ini : bool)
(notation : string) :
DC.sec_entities * unit with { title = title } =

let title = title _req_params
and rec_eq, pars = diffeqtorec eqn yofz uofk ini in

(* Not used because of following FIXME.
(* The second member of diffeqtorec's return value ist firstk,
the index from which on the recurrence holds. *)
let rec_eq = << op(1,[\$(symb:rec_eq)]) >> in
*)

(* FIXME: This has to become a test file. *)
(* FIXME: Such comments, <!-- ... --> cannot be expressed. *)
(*
if <:bool< \$(symb:rec_eq) = \$(symb:desired_result) >> then
[t: <!-- diffeqtorec passed the test for sf_id=\$(symb:sf_id)
and eqn=\$(symb:eqn). --> :t]
else
[t: <p class='warning'>Proof failed:
The return value of diffeqtorec_ddmf_wrapper<br/>\$(symb:rec)<br/>
does not match the desired result<br/>\$(symb:desired_result). :t]
endif ;
*)

(DC.section title pars, ())
```

Generated by GNU Enscript 1.6.5.90.