# $Source: /u/maple/research/lib/DEtools/diffop/src/RCS/formal_sol,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca $ # Procedures for users # formal_sol, but I think I want to write a user interface for that, so then # this procedure would fall in the next category: # Procedures for use only in the rest of the diffop package: # none # Procedures for use only in this file: # solve_semiregular, integrate_logs, lift_integral, sol_1order_eq, # lift_sol1 ############################################### # Formal solutions of differential equations #-------------------------------- ############################################### diffop_formalsol # This part of diffop computes solutions of f=xDF^n+laur*xDF^(n-1)+... # in algclos( k((x)) )[log(x),exp( .. )] # All solutions are given up to conjugation over k((x)) # The result is a list which contains a list of solutions for every exponential # part e. Such a list for 1 exp part looks like: # [ [a1,a2,..] , e] where a.i in k((x))[e,log(x)] # The corresponding solutions are: a.i*ExpInt(e) where ExpInt(e)=exp(int(e/x,x)) # If any ramifications (then x is not really x but a root of a*t^n=x) have been # used then the list is longer and looks like: # [ [a1,a2,..] , e,ext,ram,`alg over k((x))` ] where ram=a*x^n # The same form is used when algebraic extensions were needed. macro( x=`DEtools/diffop/x`, DF=`DEtools/diffop/DF`, xDF=`DEtools/diffop/xDF`, g_ext=`DEtools/diffop/g_ext`, set_laurents=`DEtools/diffop/set_laurents`, modm=`DEtools/diffop/modm`, g_expand=`DEtools/diffop/g_expand`, accuracy_laurent=`DEtools/diffop/accuracy_laurent`, value_laurent=`DEtools/diffop/value_laurent`, new_laurent=`DEtools/diffop/new_laurent`, nmterms_laurent=`DEtools/diffop/nmterms_laurent`, eval_laurent=`DEtools/diffop/eval_laurent`, formal_sol=`DEtools/diffop/formal_sol`, solve_semiregular=`DEtools/diffop/solve_semiregular`, integrate_logs=`DEtools/diffop/integrate_logs`, lift_integral=`DEtools/diffop/lift_integral`, sol_1order_eq=`DEtools/diffop/sol_1order_eq`, lift_sol1=`DEtools/diffop/lift_sol1` # test_result=`DEtools/diffop/test_result` ): # :) macro( l_p=`DEtools/diffop/l_p`, factor_op=`DEtools/diffop/factor_op` ): # Adapted to Maple 5.4 syntax for [] # Input: an operator f of order >= 1 # Output: the formal solutions formal_sol:=proc(f,ext) global x,DF; local v,i,res,c,lp; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if nargs=1 then return procname(f,g_ext(f)) fi; if has(f,DF) then return procname(l_p(f,0),ext) fi; v:=factor_op(f,'semireg',ext); res:=NULL; for i in v do c:=solve_semiregular(factor_op(i[1],`split over k((x))`,i[2]),i[2]); lp:=[c,i[4]]; if not (i[2]=ext and i[3]=x) then lp:=[op(lp),op(i[2..3]),`alg over k((x))`] fi; # test_result(formal_sol,f,ext,op(lp)[1..2],i); res:=res,lp od; [res] end: # Input: a list of factors [xDF+r1,...xDF+rn] where # the 0'th coefficients of r.i are in Z. # Output: the solutions solve_semiregular:=proc(v,ext) global x,xDF; local s,i; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if v=[] then return [] fi; s:=procname([seq(xDF+eval_laurent(v[i]-v[nops(v)] ,ext),i=1..nops(v)-1)],ext); s:=[op(map(integrate_logs,s,ext)),1]; [seq(eval_laurent(expand(i*sol_1order_eq(coeff(v[nops(v)],xDF,0) ,ext)),ext,log(x)),i=s)] end: # used by solve_semiregular # determine g such that x*d/dx (g)=ff where ff in k((x))[log(x)] integrate_logs:=proc(ff,ext) global x; local i,d,l,f; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; d:=degree(ff,log(x)); if d=0 then if has(ff,set_laurents) then f:=eval_laurent(ff,ext); # l becomes the integral of f/x: nmterms_laurent(f,1,0)*log(x)+new_laurent(1,expand(int( nmterms_laurent(f,0)/x,x)),[lift_integral,f,ext]) else g_expand(int(ff/x,x),ext) fi else # Bugfix (oct 94) l:=procname(coeff(ff,log(x),d),ext); l:=coeff(l,log(x),0)+log(x)*coeff(l,log(x),1)/(d+1); expand(l*log(x)^d+procname(convert([seq(coeff(ff,log(x),i) *log(x)^i,i=0..d-1)],`+`)-d*log(x)^(d-1)*coeff(l,log(x),0),ext)) fi end: # Used by integrate_logs # lifts LL which is int(aa/x , x) lift_integral:=proc(LL,aa,ext,acc) global x; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; value_laurent[LL]+expand(int( nmterms_laurent(aa,acc,accuracy_laurent[LL])/x,x)) end: # used by solve_semiregular # Input: a Laurent series a in Z + x k[[x]] # Output: a solution L of a*L + x*d/dx L = 0 sol_1order_eq:=proc(a,ext) global x; local n; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; n:=-nmterms_laurent(a,1,0); # n := coefficient of x^0 in a if not type(n,integer) then error "should not happen" fi; new_laurent(n+1,x^n,[lift_sol1,a,n,ext]) end: # used by sol_1order_eq lift_sol1:=proc(LL,aa,n,ext,acc) global x; local L,a,t,i; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; L:=value_laurent[LL]; a:=nmterms_laurent(aa,acc-n); for t from accuracy_laurent[LL] to acc-1 do # divide the coefficient of x^t in a*L by n-t: L:=L+modm(g_expand(convert([seq(coeff(L,x,i)*coeff(a,x,t-i) ,i=ldegree(L,x)..degree(L,x))],`+`)/(n-t)*x^t,ext)) od; L end: #savelib('formal_sol','solve_semiregular','integrate_logs','lift_integral',\ 'sol_1order_eq','lift_sol1'):