# $Source: /u/maple/research/lib/DEtools/src/RCS/expsols,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca $ # $Notify: glabahn@daisy.uwaterloo.ca $ # # DEtools[expsols] computes for the linear ordinary differential # equation # # An*y^(n) + .. A1*y' + A0*y = g (1) # # a basis for the space of exponential solutions y of (1) i.e. solutions of the # type # # y(x) = exp(int(f(x),x)) # # where the Ai's are polynomials with coefficients in some constant field K, # and where f is a rational function in x over the algebraic closure of K. # The exponential solution are used for an order reduction. # # # CALLING SEQUENCE : # # DEtools[expsols](A,g,x); # # input : A - list of the coefficient polynomials in x where A[i] is # the coefficient of the (i-1)th derivative, i.e. # A = [A0,A1,..,An] with the Ai's from (1) # g - nonhomogeneous part (arbitrary expression in x) # x - variable # # output : homogeneous case (g=0) : # [ basis functions ] # nonhomogeneous case : # [ [ basis functions ], particular solution ] # # If not all basis functions can be found then the last element # of the basis function list contains an unevaluated DEsol- # expression as well as the particular solution if g <> 0. # # NOTE: the reduction of order in this file will get obsolete and then # get removed. Use DEtools[reduceOrder] instead. macro( dsolve_int = `dsolve/int`, polylinearODE = `dsolve/diffeq/polylinearODE`, reduceorder = `DEtools/expsols_reduceorder`, simplifybasis = `dsolve/diffeq/linearODE_simplifybasis`, completesolutions = `dsolve/diffeq/linearODE_completesolutions`, ispolylinearODE = `dsolve/diffeq/ispolylinearODE`, expsols = `DEtools/expsols`, index_them = `DEtools/diffop/index_them`, split_them = `DEtools/diffop/split_them`, degree_Z = `DEtools/diffop/degree_Z`, subs_Z = `DEtools/diffop/subs_Z` ): # # expsols(A,g,x): exponential type solutions of the LODE given by A # # An*y^(n) + .. A1*y' + A0*y = g -> A = [A0,A1,..,An] # # input : A - list of coefficient polynomials in x where A[i] is the # coefficient of the (i-1)th derivative # g - right hand side of the l.o.d.e. (arbitrary expression in x) # x - variable # # output : homogeneous case (g=0) : # [ basis functions ] # nonhomogeneous case : # [ [ basis functions ], particular solution ] # expsols:= proc(A,g,x) local sols,order,answer, AA, gg; option `Copyright (c) 1997 by the University of Waterloo. All rights reserved.`; `ODEtools/userinfo`(3,'dsolve',`trying exponential solutions`); if not type(A,list) then answer := DEtools[convertAlg](A,g); if not answer = FAIL and ispolylinearODE(answer[1],answer[2],op(g),'AA','gg')then return procname(AA,gg,op(g)) else return [] fi; elif not type(A,list(polynom(anything,x))) then if ispolylinearODE(A,g,x,'AA','gg')then return procname(AA,gg,x) else error "Differential equation needs to have rational function coeffs"; fi; elif g = 0 and not has( A, x) then _EnvExponentialForm := true: return DEtools[constcoeffsols](A,x) elif g = 0 then _EnvExponentialForm := true: answer := DEtools[eulersols](A, x) ; if answer <> [] then return answer fi; fi; # should we be doing reduction of order or ONLY finding exponential # solutions ? ( false = only, true = reduce) if nargs = 4 then answer := args[4]; else answer := false; fi; sols:=`DEtools/Expsols`(A,g,x); order:= nops(A)-1; if g<>0 then gg:=`if`(nops(sols)>1,sols[2],NULL); sols:=sols[1] fi; if nops(sols) = order then `ODEtools/userinfo`(2,'dsolve',`exponential solutions successful`) elif sols <> [] then `ODEtools/userinfo`(2,'dsolve',`expon. solutions partially successful. Result(s) =`, sols) elif g <> 0 and gg <> NULL then `ODEtools/userinfo`(2,'dsolve',`particular exponential solution found`, gg) fi; sols:= simplifybasis(sols,x); if not answer or (nops(sols)=order and gg<>NULL) then # return only exponential solutions if g <> 0 then [sols,gg] else sols fi; elif g=0 and sols <> [] and nops(sols) < order then sols:= reduceorder(A,0,sols,x); simplifybasis(sols,x) elif g<>0 and sols <> [] and nops(sols) < order then sols:= reduceorder(A,0,sols,x); if g <> 0 then if gg<>NULL then [sols,gg] else [sols,`DEtools/ratsols/particular`(A,g,x)] fi else sols fi else completesolutions(A,g,sols,x) fi end: # # solvericcat(A,maxsols,x) : rational solutions of the assoc. Riccati equation # # input : A - coefficient polynomials # maxsols - maximum number of solutions to find # x - variable # # output : list of all the rational solutions of the associated Riccati eq. # # solvericcati:= proc(A,x) # local ans,i; # option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; # ans:=`DEtools/Expsols`(A,0,x,`use Int`); # # remove the exp(Int( ..)) # ans:=[seq(evala(Normal(diff(i,x)/i)),i=ans)]; # if has([args],`no indexed RootOfs`) then # map(split_them,ans,A) # else # map(index_them,ans,A) # fi # end: # Compute all conjugates of the RootOf's in the splitting field. split_them:=proc(i,A) local j,v,w,T; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; v:=indets(i,RootOf) minus indets(A,RootOf); if v={} then return i fi; while indets(map(op,v),RootOf)<>{} do v:=indets(map(op,v),RootOf) od; v:=v[1]; w:=[seq(RootOf(j[1],T),j=polytools[splits](subs_Z(op(v),T),T)[2])]; seq(procname(subs(v=j,i),{A,j}),j=w) end: # Replace the RootOf's that don't occur in A by indexed RootOf's. # Currently not used. #index_them:=proc(i,A) # local j,v; # option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; # v:=indets(i,RootOf) minus indets(A,RootOf); # v:={seq(`if`(has(j,index),NULL,j),j=v)}; # if v={} then # return i # fi; # v:=v minus indets(map(op,v),RootOf); # v:=v[1]; # seq(procname(subs(v=RootOf(op(v),index=j),i),A),j=1..degree_Z(op(1,v))) #end: # Currently not used. #degree_Z:=proc(i) # local j; # option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; # degree(subs_Z(i,j),j) #end: subs_Z:=proc(i,x) global _Z; local j,v; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; v:=indets(i,RootOf); if v={} then subs(_Z=x,i) else subs(j=v[1],procname(subs(v[1]=j,i),x)) fi end: # # reduceorder(A,g,partsols,x) # # input : A - coefficient polynomials # g - right hand side (arbitrary expression in x) # partsols - nonempty list of partial (expon. type) solutions # x - variable # # output : homogeneous case (g=0) : # [ basis functions ] # nonhomogeneous case : # [ [ basis functions ], particular solution ] # # NOTE : if the partial solutions are not of exponential type, i.e. not of the # type exp(int(f,x)) where f is a rational function, then the # coefficients of the reduced equations are in general no polynomials # and the (recursive) call of polylinearODE is not valid. # reduceorder:= proc(A,g,partsols,x) local sols,y,dy,D,rest,order,comgcd,comlcm,comfac,i,j,nrofpartsols; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; order:= nops(A)-1; `ODEtools/userinfo`(3,'dsolve',`reduction to order`,order-1); y:= partsols[1]; if order = 2 then sols:= [y,y*int(simplify(expand(exp(-int(A[2]/A[3],x)))/y^2,exp),x)]; `ODEtools/userinfo`(3,'dsolve',`back in order`,order); if g = 0 then return sols else rest := completesolutions(A,g,sols,x); return [sols,rest[2]] fi; fi; dy[0]:= y; for i from 1 to order-1 do dy[i]:= diff(dy[i-1],x) od; D:= [ seq( normal( add(A[i+1] * i!/(j+1)!/(i-j-1)! * dy[i-j-1], i=j+1..order) ), j=0..order-1) ]; comlcm:= frontend(lcm,[seq(denom(D[i]),i=1..order)]); comgcd:= numer(D[1]); for i from 2 to order while comgcd <> 1 do comgcd:= frontend(gcd,[comgcd,numer(D[i])]) od; comfac:= normal(comlcm/comgcd); D:= [seq(normal(comfac*D[i],expanded),i=1..order)]; nrofpartsols:= nops(partsols); if nrofpartsols > 1 then rest:= [ seq(simplify(diff(partsols[i]/y,x)), i=2..nrofpartsols) ]; sols:= reduceorder(D,normal(comfac*g),rest,x); else sols:= polylinearODE(D,normal(comfac*g),x) fi; `ODEtools/userinfo`(3,'dsolve',`back in order`,order); if g = 0 then [op(partsols),seq(y*dsolve_int(sols[i],x), i=nrofpartsols..nops(sols)) ] else [ [op(partsols), seq(y*dsolve_int(sols[1][i],x), i=nrofpartsols..nops(sols[1])) ], y*dsolve_int(sols[2],x) ] fi end: #savelib('expsols'): #savelib('reduceorder'): macro( DF=`DEtools/diffop/DF`, x=`DEtools/diffop/x`, expsols=`DEtools/diffop/expsols`, to_diffop=`DEtools/diffop/to_diffop`, to_eqn=`DEtools/diffop/to_eqn` ): # Input: a differential equation f. # Output: exponential solutions `DEtools/Expsols`:=proc(f,dvar) global DF,x; local L,i,v; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; `DEtools/diffop/l_p`(); if not has({args},`use Int`) then L:=procname(args,`use Int`); `ODEtools/userinfo`(3,'dsolve',nops(L),`exponential solutions found`); if not assigned(_Env_odsolve_int) then _Env_odsolve_int:=int fi; if nops(L)>0 and _Env_odsolve_int<>Int then userinfo(3,'diffop', `Simplifying exponential integral with`,_Env_odsolve_int); simplify(expand(subs(Int=_Env_odsolve_int,L))) else L fi elif indets([args],{radical,nonreal})<>{} then L:=[args]; v := radfield(indets(L,{radical,nonreal})); eval(subs(v[2],procname(op(eval(subs(v[1],L)))))) elif not has({args},`no conjugates`) then i:=procname(args,`no conjugates`); if i<>[] and type(i[1],list) then [map(split_them,i[1],{args}),op(i[2..nops(i)])] else map(split_them,i,{args}) fi elif type(f,list) then # subs( options = NULL ) so the options remember works. L:=subs(args[3]=x,`use Int`=NULL,`no conjugates`=NULL, [add(f[i]*DF^(i-1),i=1..nops(f)), args[2],args[4..nargs]]); subs(x=args[3],expsols(op(L))) elif type(dvar,function) then procname(op(DEtools[convertAlg](args[1..2])),op(dvar),args[3..nargs]) else error fi end: #savelib('`DEtools/Expsols`',\ 'split_them','subs_Z'):