# $Source: /u/maple/research/lib/DEtools/src/RCS/expsols,v $ # $Notify: hoeij@sci.kun.nl $ # # 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 = readlib(`dsolve/int`), polylinearODE = `dsolve/diffeq/polylinearODE`, reduceorder = `DEtools/expsols_reduceorder`, simplifybasis = readlib('`dsolve/diffeq/linearODE_simplifybasis`'), completesolutions=readlib('`dsolve/diffeq/linearODE_completesolutions`'), ispolylinearODE = readlib('`dsolve/diffeq/linearODE_ispolylinearODE`'), expsols = `DEtools/expsols` ); macro( 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.`; 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 RETURN( DEtools[constcoeffsols](A,x) ) elif g = 0 then 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:=readlib(`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 userinfo(2,'dsolve',`exponential solutions successful`) elif sols <> [] then userinfo(2,'dsolve',`expon. solutions partially successful. Result(s) =`, sols) fi; sols:= simplifybasis(sols,x); if not answer or (nops(sols)=order and gg<>NULL) then # return only exponential solutions if g <> 0 and sols <> [] 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 readlib(`DEtools/ratsols`): [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:=readlib(`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=readlib(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; 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)]; 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( convert([ seq(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; 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',\ 'reduceorder',\ '`DEtools/expsols.m`'): #savelib('reduceorder', '`DEtools/expsols_reduceorder.m`'): macro( DF=`DEtools/diffop/DF`, x=`DEtools/diffop/x`, expsols=`DEtools/diffop/expsols`, readlibs=readlib(`DEtools/diffop/readlibs`), to_diffop=readlib(`DEtools/diffop/to_diffop`), to_eqn=readlib(`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.`; # load all files used by the expsols code: readlibs(expsols); if not has({args},`use Int`) then L:=procname(args,`use Int`); userinfo(3,'diffop',`Expanding the integral symbols`); simplify(expand(subs(Int=int,L))) elif indets([args],radical)<>{} then L:=[args]; L:=procname(op(subs(readlib(`radnormal/radfield2`)(L,'v'),L))); subs(eval(v),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 L:=subs(args[3]=x,[convert([seq(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() # L:=traperror(to_diffop(args)); # if L=lasterror then # ERROR(L) # fi; # [seq(to_eqn(i,args[2..nargs]),i=expsols(L))] fi end: #savelib('`DEtools/Expsols`',\ 'split_them','subs_Z','`DEtools/Expsols.m`'):