# $Source: /u/maple/research/lib/DEtools/src/RCS/kovacicsols,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca # # #--> kovacicsols: Kovacic's algorithm for solving second order o.d.e.'s # # Calling sequence: kovacicsols(A,var) # # Purpose: Solve a second order linear homoegeneous differential # equation of the form # # A[3] * y'' + A[2] * y' + A[1] * y = 0 # # Input: A - list of coefficient polynomials in var # # Output: list with two Liouvillian solution, if such solutions # exist, empty list otherwise. # # interface adapted, S.Schwendimann April/93 # adapted for use in DEtools, G. Labahn Sept/96 # # Completely rewritten, Mark van Hoeij, May/2000. macro( REDUC = 1, # Reducible without algebraic extension, Cases 1 and 2a. R_IMPR = 2, # Reducible or imprimitive. Cases 2b and 3. D2 = [4,{1,2}], # D2_SL2 = D4 with 2nd degree Riccati defined over # algebraic extension of degree 3. Case 4. A4 = [4,{1,2,3}], # A4_SL2 = tetrahedral group, with degree 4 Riccati # defined over base field. Case 5. A4alg = [6,{1,2,3}], # Case 6, a tetrahedral group with degree 4 # Riccati defined over quadratic extension # computed by subfields(S2,2,ext,x). S4 = [6,{1,2,3,4}], # Case 7, octahedral group S4_SL2. A5 = [12,{1,2,3,5}] # A5_SL2 = icosahedral group, case 8. ): # What the code currently does in various cases: # ---------------------------------------------- # # case 1) reducible and the group is not reductive # # L = L1*L2 and L2 is the only right-hand factor, so there is only 1 # exponential solution, exp(int(...)). If that solution looks # complicated, then it uses Int to find the other solution, otherwise # it uses int. # # case 2a) reducible without algebraic extensions, group is reductive # # The code computes two exponential solutions, it uses exp(int(...)). # # case 2b) L is reducible over an algebraic extension, group is reductive # # The code computes two exponential solutions, it uses exp(int(...)). # # case 3) group is imprimitive and the Riccati polynomial of degree two is # defined over the base field. # # This boils down to finding a first order factor over the base field # of the 2nd symmetric power. Once you have that, then after a # substitution DF -> DF + something the 2nd symmetric power has a # right-hand factor DF. When the 2nd symmetric power has DF as right # factor then it can be solved by plus_min_expsols. # # case 4) group is imprimitive and the Riccati polynomial of degree # two is NOT defined over the base field. # # This means that the group is D2 (notation from Ulmer/Weil paper). # Initially the code finds a Riccati polynomial of degree 4 very # quickly. The method in the Ulmer/Weil paper, implemented on: # http://www-lmc.imag.fr/lmc-cf/Manuel.Bronstein/kovacic_demo.html # by Bronstein, stops at this point and returns the degree 4 # polynomial. However, a quadratic Riccati polynomial is definitely # nicer than one of degree 4 because you can express the answer with # a nice radical instead of a RootOf that users are not so fond of. # Therefore, this code does considerably more effort (read spends # more CPU time) on this case. It calculates the cubic extension # over which the degree 2 Riccati polynomial is defined, and then # runs the code again, restricting the search to Riccati polys of # degree 2). # This default can be altered by setting: _EnvkovacicsolsD2:=false # With that setting the code will behave like Ulmer/Weil on case 4. # But case 3 will still be different, the Ulmer/Weil method does # not search for exp sols of the 2nd symmetric power, hence it may # not find the degree 2 Ricatti poly in case 3, whereas this code # does. # # case 5) group is A4 and the Riccati polynomial of degree 4 is defined # over the base field. # # The Ulmer/Weil method will not necessarily find this degree 4 Riccati, # because it boils down to searching exponential solutions. It will usually # find the one of degree 6. When I tried Bronsteins code it also returned # a poly of degree 6, so it follows Ulmer/Weil's paper. # This code does more effort, if there is a Riccati of degree 4 over # the base field, it finds it. # # case 6) group is A4 and the Riccati polynomial of degree 4 is NOT defined # over the base field. # # If _EnvkovacicsolsA4<>false (the default) then it will compute the # quadratic field extension over which a degree 4 Riccati can be found, # and runs the code again, restricting the search to degree 4 Riccati. # Going from degree 6 to degree 4 is not always worth introducing a # quadratic extension of the constants, so this default can be changed # by setting: _EnvkovacicsolsA4:=false; then this code follows # Ulmer/Weil, it computes a Riccati polynomial of degree 6. # # case 7) S4 # # The code returns a Riccati polynomial of degree 6, just # like the original Kovacic algorithm. Ulmer/Weil uses degree 8 here. # # case 8) A5 # # The code returns a Riccati polynomial of degree 12, like Ulmer/Weil and # like the original Kovacic algorithm. # # # # Summary: # -------- # # Essentially this code computes the same as the original Kovacic # algorithm. But some defaults can be changed: # If _EnvkovacicsolsA4=false then in case 6 it follows Ulmer/Weil. # If _EnvkovacicsolsD2=false then in case 4 it follows Ulmer/Weil. # # This code should not bog down when the singularities are # not rational. This code handles those field extensions much more # carefully than the original Kovacic algorithm and so it should perform # much better in those cases. # To make the outputs like nice for the imprimitive groups the # algorithm chooses carefully between int and Int. # # # What remains to be done: # ------------------------ # # 1) What still needs to be done for order 2 is to simplify the # exp(Int(RootOf( Riccati polynomial ))) for the cases 4)-8). # This not only gets rid of the exp and the Int, but it also leads # to a minimum polynomial with smaller coefficients. So this is # worthwhile. There are several different methods for doing this. # # Done so far: haven't put this in yet, but instead for A4 and S4 # I try to get solutions in terms of radicals when it # is not too difficult to get them. # # 2) Do something about order 3. At least do the case where # the Riccati poly has degree 3 (imprimitive group) and # is defined over the base field. That shouldn't be too hard. # # Done so far: there is some experimental code by J.A. Weil and # me that we never finished. Should do that some time. macro( order2_gen_exp =`DEtools/kovacicsols/gen_exp`, order2 =`DEtools/kovacicsols/order2`, symproduct_order2=`DEtools/kovacicsols/sympr`, sympower_order2 =`DEtools/kovacicsols/sympow`, sympowsol =`DEtools/kovacicsols/sympowsol`, seri =`DEtools/kovacicsols/seri`, find_answer =`DEtools/kovacicsols/find_answer`, which_group =`DEtools/kovacicsols/which_group`, construct_riccati=`DEtools/kovacicsols/constr_r`, e_int =`DEtools/kovacicsols/e_int`, combs =`DEtools/kovacicsols/combs`, g_roots =`DEtools/kovacicsols/roots`, g_factors =`DEtools/kovacicsols/factors`, plus_min_expsols =`DEtools/kovacicsols/plus_min_expsols`, combi =`DEtools/kovacicsols/combi`, simpl_prod =`DEtools/kovacicsols/simpl_prod`, transfo =`DEtools/kovacicsols/transfo`, convert_ln =`DEtools/kovacicsols/convert_ln`, int_ =`DEtools/kovacicsols/int`, real_answer =`DEtools/kovacicsols/real_answer`, EXT = {'nonreal', 'radical', 'algext'}, ispolylinearODE = `dsolve/diffeq/ispolylinearODE`, userinfo = `ODEtools/userinfo`, sort_poly =`polytools/sort_poly`, subfields =`evala/Subfields` ): `DEtools/kovacicsols` := proc(A,var) local answer, AA, gg,i; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: M. van Hoeij`; # Standard conversion equation -> list: if not type(A,list) then answer := DEtools['convertAlg'](A,var); if not answer = FAIL and ispolylinearODE(answer[1],answer[2],op(var),'AA','gg') then if gg = 0 then return procname(AA,op(var),args[3..nargs]); else error "Input should be a homogeneous ODE" fi; else error "Input should be linear ODE with polynomial coefficients" fi; fi; if nops(A) <> 3 then error "Only second order is implemented" fi; gg:=add(A[i]*AA^i,i=1..nops(A)); if indets(map(denom,A),'name')<>{} or has(content(gg,AA),var) then gg:=evala(Primpart(gg,AA)); return procname([seq(coeff(gg,AA,i),i=1..nops(A))],args[2..-1]) fi; answer:=order2(args); if answer<>[] and _Envkovacicsols_check=true then # Set this Env variable for testing and debugging. Digits:=24; for AA in answer do gg:=A[1]*AA+A[2]*diff(AA,var)+A[3]*diff(AA,var,var); if has(gg,{'Int','int'}) and simplify(evala(numer(gg)))=0 then next fi; gg:=subs([seq(i=rand(210..290)()/100.0,i=indets(gg,'name'))],gg); gg:=evalf(gg); if abs(gg) < 0.1^(Digits/2) then next fi; error `check failed` od; userinfo(1,'dsolve',`answer verified`) fi; answer end: # A is list of polynomials # if p is a RootOf it should also appear in ext. order2_gen_exp := proc(A::list, T::name, x::name , p, ext::set) local B,i,DF,m,a0,a1,a2,a3,b0,b1,c0; Normalizer:=`if`(ext={} and not has(p,'RootOf'),normal,`evala/Normal`); if nops(A)<>3 or A[3]=0 then error "wrong input" elif p=infinity then B:=[A[1]*x^4, collect(-A[2]*x^2+2*x*A[3],x,Normalizer), A[3]]; m:=max(op(map(degree,B,x))); B:=[seq(collect(x^m*subs(x=1/x,i),x),i=B)]; procname(B,T,x,0,ext) else B:=[seq(series(subs(x=x+p,i[1]),x=0,i[2]),i=[[A[3],4],[A[2],3]])]; a0,b0 := op(map( Normalizer, map(coeff,B,x,0)) ); if a0<>0 then # `DEtools/diffop/alert`:="nonsingular"; [0,1,T=x] else a1:=Normalizer(coeff(B[1],x,1)); c0:=Normalizer(subs(x=p,A[1])); if b0=0 and c0=0 then error "input polys should not have common factor" elif a1<>0 then [0,Normalizer(1-b0/a1),T=x] else a2:=Normalizer(coeff(B[1],x,2)); b1:=Normalizer(coeff(B[2],x,1)); if a2<>0 and b0=0 then B:=a2*x^2+b1*x-a2*x+c0; i:=g_roots(B,x,ext); if i<>[] then [op(i),T=x] else [RootOf(B,x),T=x] fi else a3:=Normalizer(coeff(B[1],x,3)); if a2<>0 then [0, Normalizer(2-b1/a2+b0*a3/a2^2) - Normalizer(b0/a2/T),T=x] elif a3<>0 and b0=0 and b1=0 then [1/T+Normalizer(3/4-coeff(B[2],x,2)/2/a3), Normalizer(-c0/a3)*T^2 = x] else # We've captured the most common cases, handle the more # difficult cases with the complete algorithm: B:=DEtools['gen_exp'](A[3]*DF^2+A[2]*DF+A[1],[DF,x],T,x=p, 'groundfield'=ext); if nops(B)=2 then [B[1][1],B[2][1],T=x] else B[1] fi fi fi fi fi fi end: g_roots:=proc(f,x::name) local v,i,B; B:=evalb(indets(f,'name')={x}); if B then v:=traperror(roots(args)) fi; if v=lasterror or not B then v:=g_factors(args); [seq(`if`(degree(i[1],x)=1,evala(Normal( -coeff(i[1],x,0)/coeff(i[1],x,1)))$i[2],NULL),i=v)] else map(i -> i[1]$i[2], v) fi end: g_factors:=proc(f,x,ext) local v; try v:=evala(Factors(f,{op(ext)})); catch : v:=evala(Factors(numer(normal(f)),{op(ext)})) end try; v:=select(has,v[2],x) end: order2:=proc(A::list, x::name, ext::set, case::list) local sing,dp,p,q,g,r,dr,i,T,pos1,pos2,pos4,pos6,posh,m,pos_set ,cases,cm,DF,S2,S3,ans,threes,twos; global `DEtools/kovacicsols/T`; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: M. van Hoeij`; if nops(A)<>3 then error "wrong arguments" elif nargs=1 then p:=indets(A,'name'); if nops(p)<>1 then x # = error fi; return procname(A,p[1]) elif nargs=2 then return procname(A,x,indets(A,EXT)) fi: Normalizer:=`if`(ext={},normal,`evala/Normal`); S2,S3,threes,twos:={}$4; cases:=[REDUC,R_IMPR,D2,A4,A4alg,S4,A5]; if nargs=4 then cases:=case fi; pos_set:={op(cases)}; sing:=sort_poly([infinity,seq(i[1],i=g_factors( A[3],x,[op(ext)]))],x); cm:=NULL; for q in sing while pos_set<>{} do if q=infinity then p:=q; dp[p]:=1 else p:=RootOf(q,x); dp[p]:=degree(q,x) fi; cm:=cm,p; g[p]:=subs(T=x,order2_gen_exp(A,T,x,p,ext union `if`(dp[p]>1,{p},{}))); if nops(g[p])=2 then # generalized exponents are defined over an extension. pos_set:=pos_set intersect {R_IMPR}; if g[p][2]=(x=x) then m:=evala(Trace(g[p][1],ext,ext)) else m:=x^2/lhs(g[p][2]); m:=add(Normalizer(coeff(g[p][1],x,2*i)*m^i)*2*x^i, i=ceil(ldegree(g[p][1],x))..0) fi; r[p]:=0; pos2[p]:={m} else m:=Normalizer(g[p][1]+g[p][2]); pos1[p]:={g[p][1],g[p][2]}; pos2[p]:={m}; r[p]:=Normalizer(g[p][1]-g[p][2]); if type(r[p],rational) then if r[p]<0 then i:=g[p][1]; r[p]:=-r[p] else i:=g[p][2] fi; if r[p]=0 then # Local solution involves logarithm pos1[p]:={i+r[p]}; pos_set:=pos_set intersect {REDUC} else dr:=denom(r[p]); if dr=1 then pos1[p]:={i}; pos2[p]:={2*i} elif dr=2 then twos:=twos union {[p,dp[p]]}; pos2[p]:={m,2*i} elif dr=3 then threes:=threes union {[p,dp[p]]}; if irem(dp[p],2)=0 then # ext of degree 2 for A4 S2:=S2 union {q} else pos_set:=pos_set minus {A4alg} fi fi; pos4[p]:={seq(2*m+(i-2)*r[p],i=0..`if`(dr=3,1,0))}; pos6[p]:={seq(3*m+(i-3)*r[p]+`if`(dr=4,r[p],0), i=0..`if`(dr=2,1,0))}; pos_set:=pos_set minus {seq(`if`(member(dr,i[2]),NULL ,i),i=pos_set minus {REDUC,R_IMPR})} fi else r[p]:=0; pos_set:=pos_set intersect {REDUC,R_IMPR} fi; if nops(pos1[p])=2 then if irem(dp[p],3)=0 then S3:=S3 union {q} else pos_set:=pos_set minus {D2} fi fi fi; posh[p]:=m/2 od; if member(D2,pos_set) then # If dr=denom(r[p]) never becomes 3 then the group can't # be A4. In fact you need at least two 3's. pos_set:=pos_set minus {A4,A4alg} fi; sing:=[cm]; cm:=0; while not type(cm,list) do cm:=[seq(`if`(member(m,pos_set), [m,combs(ext,x,[seq([p,dp[p], `if`(m=REDUC,pos1[p], `if`(m=R_IMPR,pos2[p], `if`(m=A4,pos4[p], `if`(m=S4,pos6[p], {m[1]*posh[p]-trunc(m[1]*r[p]/2)})))) ],p=sing)])],NULL),m=cases)]; m:=add(op(1,i[1])*nops(i[2]),i=cm); if m>40 and pos_set<>{REDUC} then for p in sing do if type(r[p],posint) and r[p]<10 and DEtools['formal_sol'](add(A[i+1]*DF^i,i=0..2) ,[DF,x],`has logarithm?`,x=p) then # there are logarithms. pos1[p]:={op(pos1[p])+r[p]}; pos_set:=pos_set intersect {REDUC}; cm:=0; break fi od fi; od; ans:=[]; for m in cm while ans=[] do for r in m[2] do g:=[`if`(r[1][3]=0,0,add(-coeff(r[1][3],x,i)*x^(-i-1), i=ldegree(r[1][3],x)..-1)), seq(subs(x=x-r[i][1],collect(r[i][3]/x,x)),i=2..nops(r)-1)]; g:=Normalizer(add(`if`(r[i][2]=1,Normalizer(g[i]),evala( Trace(normal(g[i]),ext,ext))),i=1..nops(r)-1)/op(1,m[1])); dr:=symproduct_order2(A,g,x); if op(1,m[1])<3 then T:=sympower_order2(dr,op(1,m[1]),x); i:=DEtools['polysols'](T,0,x) else i:=sympowsol(dr,r[-1],op(1,m[1]),x) fi; if i<>[] then userinfo(1,'dsolve',`A Liouvillian solution exists`); if m[1]=S4 then # Set of points where the local monodromy is an # odd permutation. twos:={seq(`if`(type(Normalizer(6*posh[i[1]]-i[3]) ,integer),NULL,i[1..2]),i=r[1..-2])} fi; g:=A,x,ext,i,g,dr,m[1],S2,S3,threes,twos; i:=traperror(find_answer(g)); if i=lasterror then if op(1,m[1]) < 3 then error i fi; i:=1/rand(2..7)(); i:=sympowsol(dr,r[-1],op(1,m[1]),x,i); if i=[] then next fi; g:=subsop(4=i,[g]); i:=find_answer(op(g)) fi; ans:=[op(ans),op(i)]; if nops(ans)=2 then if m[1]=REDUC and ans[2] = eval(subs( exp=proc(a) exp(-a) end,ans[1])) then ans:=eval([subs(exp=sinh,ans[1]) ,subs(exp=cosh,ans[1])]) fi; if member(m[1],[REDUC,R_IMPR]) and not has(ans,'RootOf') and not has(ans,'Int') and subs(I=-I,ans[1])=ans[2] and not has( indets(ans,'radical'),x) then # translate (a+b*i)*exp(c+d*i) # into sin and cos: ans:=real_answer(ans) fi; return ans elif m[1]<>REDUC then error "failed" fi fi; od od; if ans=[] then userinfo(1,'dsolve',`No Liouvillian solutions exists`); [] else userinfo(1,'dsolve',`Group is reducible, not completely reducible`); i:=e_int(`DEtools/kovacicsols/T`,x); if has(convert(i,ln),ln) or (indets(i,EXT)<>{} and not (type(i,ratpoly(anything,x)) and degree(denom(i),x)<2)) then i:=Int(i,x) else i:=int_(i,x) fi; [ans[1], combi(combine(ans[1]*i,exp),x)] fi end: # formula for symmetric_product(A[3](x)*DF^2+A[2](x)*DF+A[1](x),DF+g(x),[DF,x]) symproduct_order2:=proc(A,g,x) map(Normalizer,[g^2+diff(g,x)+A[1]/A[3]+A[2]/A[3]*g,2*g+A[2]/A[3],1]) end: sympower_order2:=proc(A,m,x) local L0,L1,i,DF; L0,L1:=1,DF; for i from 1 to m do L0,L1:=L1,collect(L1*(DF+i*A[2])+diff(L1,x) +i*(m-(i-1))*A[1]*L0,DF,Normalizer) od; L1:=primpart(L1,DF); # L1:=evala(Primpart(L1,DF)); [seq(coeff(L1,DF,i),i=0..m+1)] end: sympowsol:=proc(B,S,m,x) local A,L0,L1,i,acc,c; if type(S,nonnegint) then A:=procname(B,add(c[i]*x^i,i=0..S) ,m,x); L0:=indets(A) intersect {seq(c[i],i=0..S)}; L1:=seq(i=0,i=L0); A:=[seq(collect(subs(i=1,L1,A),x,Normalizer),i=L0)]; else acc:=degree(S,x)+3+m; # Note: this procedure is a heuristic in the sense that it # can return too many "solutions". But a check is not necessary # because correctness is checked in construct_riccati, and if # it fails then it has a mechanism to repair it. A:=[B[1]/B[3],B[2]/B[3]]; if nargs>4 or member(0,map(Normalizer,subs(x=0,map(denom,A)))) then i:=1; if nargs>4 and type(args[5],rational) then i:=args[5] fi; return subs(x=x+i,procname(subs(x=x-i,B),S,m,x)) fi; A:=map(seri,A,x,acc); L0,L1:=S,diff(S,x); for i from 1 to m do L0,L1:=L1,seri(L1*i*A[2]+diff(L1,x) +i*(m-(i-1))*A[1]*L0,x,acc-i) od; i:=`solve/linear`({coeffs(L1,x)},indets(S) minus indets([B,x])); if i=NULL then error "failed" else subs(i,S) fi fi end: seri:=proc(f,x,acc) collect(convert(series(f,x=0,acc),polynom),x,Normalizer) end: find_answer:=proc(A,x,ext,psol,g,As,m,S2,S3,threes,twos) local S,S1,L,T,DF,i; global `DEtools/kovacicsols/T`; userinfo(2,'dsolve',which_group(m)); S:=e_int(g,x); S1:=combi(S*psol[1]^(1/op(1,m)),x); L:=symproduct_order2(As,diff(psol[1],x)/psol[1]/op(1,m),x); if m=REDUC then if L[1]<>0 then error "failed" elif nops(psol)=2 then [S1,S*psol[2]] else # Set this variable to help find the second solution by # reduction of order if it is not found in another way `DEtools/kovacicsols/T`:=-L[2]; [S1] fi elif m=R_IMPR then # reducible over an algebraic number or algebraic function # extension of degree 2. S:=plus_min_expsols(L,x,S1); if S=[] then error "failed" fi; T:=combine([S1*S[1],S1*S[2]],exp); S:=map(combi,T,x); while length(S)2 or S3={} then error "failed" elif _EnvkovacicsolsD2<>false then for S in subfields(S3,3,ext,x) do userinfo(2,'dsolve',`searching over the cubic extension`,S); # Now calculate the degree 2 Riccati and solve: T:=order2(A,x,ext union {S},[R_IMPR]); if T<>[] then if not has(op(1,S),'RootOf') and nops([coeffs(op(1,S),_Z)])=2 then T:=subs(S=convert(S,radical),T) fi; return T fi od; error "failed" fi elif m=A4alg and _EnvkovacicsolsA4<>false then for S in subfields(map(primpart,S2,x),2,ext,x) do # check for A4 with an algebraic extension of degree 2 userinfo(2,'dsolve',`checking A4 over the extension`,S); T:=order2(A,x,ext union {S},[A4]); if T<>[] then if not has(op(1,S),'RootOf') then T:=subs(S=convert(S,radical),T) fi; return T fi od; error "failed" elif m=A4 then S:=add(i[2],i=threes); if S<2 then error "failed" elif S=2 then # Then we can reduce A4 to D2 by a cubic extension # of genus 0. if nops(threes)=2 then # Can do it over the rational numbers S:=[threes[1][1],threes[2][1]]; if S[2]=infinity then S:=[S[2],S[1]] fi; L:=[seq(`if`(T[2]=1,T[1],NULL),T=twos)]; # If you have a rational point with denom 2 then you're # sure we don't need "D2" (i.e. alg ext of deg 3 for # the imprimitive case), otherwise try both: S1:=x,ext,[R_IMPR,`if`(L=[],D2,NULL)]; DF:=proc(v) local x; sort(v, (i,j) -> `polytools/shorter`(i,j,x))[1] end: if S[1]=infinity then S:=S[2]; if L=[] then L:=1 else L:=DF([seq(Normalizer(1/(S-T)),T=L)]) fi; T:=order2(transfo(A,x,x^3/L+S,x),S1); return subs(x=(x*L-S*L)^(1/3),T) else if S[1]=0 then S:=[S[2],S[1]] fi; if L=[] or member(infinity,L) then L:=1 else L:=DF([seq(Normalizer((T-S[2])/(T-S[1])),T=L)]) fi; T:=order2(transfo(A,x, (x^3*S[2]-S[1]*L)/(x^3-L),x),S1); L:=[numer(L),denom(L)]; return subs(x= ((L[1]*x-L[1]*S[1])/(x*L[2]-S[2]*L[2]))^(1/3),T) fi # else # can do it over a quadratic extension. It is not # necessary to implement this because A4alg can take # care of introducing that quadratic extension. fi fi elif m=S4 then # The number of points where the permutation is odd: S:=add(S[2],S=twos); if S=0 or type(S,odd) then # If a product of permutations is 1 then you # must have an even number of odd permutations error "failed" elif S=2 then # Then we can reduce S4 to A4 by a quadratic extension # of genus 0. if nops(twos)=2 then # Can do it over the rational numbers if nops(threes)=1 and threes[1][2]=1 then # Can afterwards also reduce A4 to D2 S:=[twos[1][1],twos[2][1]]; if S[2]=infinity then S:=[S[2],S[1]] fi; L:=threes[1][1]; S1:=x,ext,[A4]; if S[1]=infinity then S:=S[2]; # Make sure to get the +/- sign right: L:=Normalizer(1/(L-S)); T:=order2(transfo(A,x,x^2/L+S,x),S1); return subs(x=(x*L-S*L)^(1/2),T) else if S[1]=0 then S:=[S[2],S[1]] fi; if L=infinity then L:=1 else L:=Normalizer( (L-S[2])/(L-S[1]) ) fi; T:=order2(transfo(A,x, (x^2*S[2]-S[1]*L)/(x^2-L),x),S1); L:=[numer(L),denom(L)]; return subs(x= ((L[1]*x-L[1]*S[1])/(x*L[2]-S[2]*L[2]))^(1/2),T) fi # else # Can not reduce A4 to D2 without going to genus>0 fi # else # Need quadratic extension of the constants. # We don't handle this case but just return the # Riccati polynomial. fi fi fi; T:=construct_riccati(L,x,m[1]); combine([S1*T[1],S1*T[2]],exp) fi; end: which_group:=proc(m) if m=A4 or m=A4alg then "Tetrahedral Galois group A4_SL2." elif m=A5 then "Icosahedral group A5_SL2." elif m=D2 then "Dihedral group D4 = D2_SL2" elif m=REDUC then "Reducible group" elif m=R_IMPR then "Group is reducible or imprimitive" elif m=S4 then "Octahedral Galois group S4_SL2" else NULL fi end: construct_riccati:=proc(A,x,m) local b,i,Y; b[0],b[1]:=1,0; for i from 1 to m do b[i+1]:=Normalizer(b[i]*i*A[2]+diff(b[i],x)+i*(m-i+1)*A[1]*b[i-1]) od; if b[m+1]<>0 then error "failed" fi; i:=add((-1)^i*b[i]*Y^(m-i)/i!,i=0..m); [exp(Int(RootOf(i,Y,'index'=1),x)), exp(Int(RootOf(i,Y,'index'=2),x))] end: # Just for rational functions, no other functions. plus_min_expsols :=proc(A,x) local ax2, i, v,c,d,dd,j,jj; if nops(A) = 3 then ax2:=Normalizer(-A[1]/A[3]); # = a(x)^2 if ax2<>0 and Normalizer(diff(ax2,x)*A[3]+ax2*2*A[2])=0 then # Would give error message if you have coeffs # that are not rational functions: v:=evala(Sqrfree(ax2,x),'expanded'); c:=mul(i[1]^floor(i[2]/2),i=v[2]); d:=mul(i[1]^(i[2] mod 2),i=v[2]); if degree(d,x)>2 # or (degree(d,x)=2 and indets(A)={x}) then # Use Maple's integrator for algebraic functions, # but only when it works, i.e., when there are no # parameters. i:=RootOf(_Z^2-d*v[1]); dd:=sqrt(d*v[1],'symbolic'); # Here I'm preventing int(sqrt(x^3+...)) to avoid # elliptic functions in the output which are sure # to lead to large answers. So I integrate RootOf( # _Z^2-(x^3+..)) instead and after that, I replace # the RootOf by sqrt, but, one should then also # replace int by Int because the args of int have # changed and I want to prevent int from being # called. i:=subs(i=dd,int_(c*i,x)); [exp(i),exp(-i)] else i:=convert_ln(int_(c*sqrt(d,'symbolic'),x)); if has([i,args],{ln,exp}) then # avoid output [sinh(c*ln(x)),cosh(c*ln(x))] but # write [x^c,x^(-c)] instead. dd:=sqrt(v[1],'symbolic'); i:=select(has,`if`(type(i,`+`),[op(i)],[i]),x); i,j:=convert(select(has,i,ln),`+`)*dd, convert(remove(has,i,ln),`+`)*dd; jj:=`if`(has(j,'RootOf'),j,radnormal(j)); if length(jj){1,-1} and type(a[1],integer) then a:=max(op(a)); return procname(mul(i[1]^(i[2]/a),i=l),x,e*a,n) fi fi fi; a:=traperror(primpart(A,x)); if a=lasterror or length(a)>length(A) then a:=A fi; if not has(a,{'RootOf','Int','int'}) then i:=radnormal(a,'expanded'); if length(i) {} then b:=collect(b,i,Normalizer); fi; b:=`if`(type(b,`+`),[op(b)],[b]); c:=select(has,b,ln); d:=remove(has,b,ln); combi(expand(exp(convert(c,`+`)))*exp(convert(d,`+`)),x,x) end: # v={ [sing::point, degext::posint, pos::set], ... } combs:=proc(ext::set,x::name,v::list) local c,w; c:=proc(s) local i,j; if nargs=1 then [seq([[s[1],s[2],i]],i=s[3])] else [seq(seq([op(i),op(j)],i=procname(s)) ,j=procname(args[2..nargs]))] fi end: w:=c(op(v)); c:=proc(s,x,ext) local i; i:=Normalizer(add(`if`(i[2]=1,coeff(i[3],x,0),evala( Trace(coeff(i[3],x,0),ext union {i[1]},ext))),i=s)); if type(i,'nonposint') then [op(s),-i] else NULL fi end: map(c, w, x, ext) end: # x=old variable. # R=new variable # z=function in R # subs x=z in A # NOTE: some code-duplication here, I should probably use PDEtools[dchange] # to do this. transfo:=proc(A,x,z,R) local S,p,q,i,v,DF; v:=subs(x=z,A); _Envdiffopdomain:=[DF,R]; S, p, q:=v[1], 1, Normalizer(1/diff(z,R))*DF; for i from 2 to nops(v) do p:=DEtools['mult'](q,p); S:=S+p*v[i] od; S:=evala(Primpart(S,DF)); [seq(coeff(S,DF,i),i=0..degree(S,DF))] end: convert_ln:=proc(a) if has(_Env_kovacicsols,`no convert`) then a else convert(a,ln) fi end: int_:=proc() if _Env_odsolve_int = 'Int' then Int(args) elif _Env_odsolve_int = 'int' then subs(int=Int,_Env_odsolve_int(args)); else `PDEtools/int`(args) fi end: real_answer:=proc(S) local loc,w,W,T,n,i; loc:=proc() local a; assume(a,'real'); a end: w:=[op(indets(S[1],'name'))]; n:=nops(w); W:=[seq(loc(i),i=1..n)]; T := subs([seq(w[i]=W[i],i=1..n)],S[1]); T := subs([seq(W[i]=w[i],i=1..n)], [Re(T),Im(T)]); if has(T,'Re') or has(T,'Im') then S else map(collect,T,exp) fi end: #savelib('order2_gen_exp', 'order2', 'symproduct_order2','g_roots',\ 'g_factors', 'find_answer', 'construct_riccati',\ 'sympowsol','seri','sympower_order2','plus_min_expsols',\ 'e_int', 'combs','which_group','`DEtools/kovacicsols`',\ 'combi','simpl_prod','transfo','convert_ln','int_','real_answer'):