# $Source: /u/maple/research/lib/dsolve/diffeq/src/RCS/higherorder,v $ # # # Authors: George Labahn and Mark van Hoeij, June 2000. macro( userinfo = `ODEtools/userinfo`, higherorder = `dsolve/diffeq/higherorder`, iszero = `dsolve/diffeq/higherorder/iszero`, ROOT = `dsolve/diffeq/higherorder/root`, is_remove_lower = `dsolve/diffeq/is_remove_lower`, remove_lower = `dsolve/diffeq/remove_lower`, is_sympow_o2 = `dsolve/diffeq/higherorder/is_sympow_o2`, solve_sympow_o2 = `dsolve/diffeq/higherorder/solve_sympow_o2`, is_LCLM_nth_root = `dsolve/diffeq/higherorder/is_LCLM_nth_root`, solve_LCLM_nth_root = `dsolve/diffeq/higherorder/solve_LCLM_nth_root`, is_LCLM_with_rat_o1 = `dsolve/diffeq/higherorder/is_LCLM_with_rat_o1`, is_LCLM_o2_same_DF1 = `dsolve/diffeq/higherorder/is_LCLM_o2_same_DF1`, solve_LCLM = `dsolve/diffeq/higherorder/solve_LCLM`, is_sympr_o2_no_DF1 = `dsolve/diffeq/higherorder/is_sympr_o2_no_DF1`, solve_sympr_o2_no_DF1= `dsolve/diffeq/higherorder/solve_sympr_o2_no_DF1`, is_sympr_o2 = `dsolve/diffeq/higherorder/is_sympr_o2`, solve_sympr_o2 = `dsolve/diffeq/higherorder/solve_sympr_o2`, to_normal_form = `dsolve/diffeq/higherorder/to_normal_form`, dint = `dsolve/int`, is_MeijerG_ode = `DEtools/MeijerGsols/is_MeijerG_ode`, MeijerG_solver = `DEtools/MeijerGsols/MeijerG_solver`, symmetric_product = DEtools['symmetric_product'], symmetric_power = DEtools['symmetric_power'], EXP_INT = `DEtools/kovacicsols/e_int`, EXT = {nonreal, radical, algext} ): higherorder := proc(A,x) local ans, AA, v, ord, i, aa, bb, shift, va, vb, FT, RAT; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Authors: G. Labahn and M. van Hoeij`; # Make A into a list if it is not already in such a form. if not type(A,list) then ans := DEtools['convertAlg'](A,x); if ans = FAIL or ans[2]<>0 then return [] else return procname(ans[1],op(x),args[3..-1]) fi fi; # Make first entry = 1 and normalize entries. Normalizer:=`if`(indets([args],EXT)={},normal,`evala/Normal`); AA := [seq(Normalizer(A[i]/A[-1]),i=1..nops(A))]; ord := nops(A)-1; FT := not evalb(nargs=4 and args[3]=`second try`); # FT = first try? RAT := type(AA,list(ratpoly(anything,x))) and (FT or not args[4]); # coeffs are rational, and (if it is the second try) were not # rational in the first try. if ord<3 then # Only methods for order >=3 in this file ans := [] elif FT and is_sympow_o2(AA,x,'v') then userinfo(2,dsolve,`Equation is the `,v[1], `-th symmetric power of`,v[2]); ans := solve_sympow_o2(v) elif is_remove_lower(AA) then userinfo(2,dsolve,`The constants are solutions, reducing order`); ans := remove_lower(AA,x) elif FT and is_LCLM_nth_root(AA,x,'v') then # A second try can only be succesful if the coeffs became constant, # which (I believe) is already checked elsewhere. userinfo(2,dsolve, `Equation has first order factor y'-R*y where R^n is rational`); ans := solve_LCLM_nth_root(v,x) elif ord=4 and is_sympr_o2_no_DF1(AA,x,'v') then userinfo(2,dsolve,`Equation is the symmetric product of` ,v[1],`and`,v[2]); ans := solve_sympr_o2_no_DF1(v) elif (FT and ord=4 and is_LCLM_o2_same_DF1(AA,x,'v')) or (RAT and is_LCLM_with_rat_o1(AA,x,'v')) then userinfo(2,dsolve,`Equation is the LCLM of`,op(v[1])); ans := solve_LCLM(v) elif (FT or RAT) and ord=4 and is_sympr_o2(AA,x,'v') then userinfo(2,dsolve,`Equation is the symmetric product of`,v[1],`and` ,v[2],op(`if`(nops(v)=4,[`and`,diff(v[3],x)-v[3]*v[4]],[]))); ans := solve_sympr_o2(v,x) elif is_MeijerG_ode(AA,x,'aa','bb','va','vb','shift') then userinfo(2,dsolve,`Equation is of hypergeometric or Meijer G type`); ans := MeijerG_solver(AA,x,aa,bb,va,vb,shift); elif FT and A[ord] <> 0 then ans := to_normal_form(AA,x); else ans := [] fi; ans; end: #savelib('higherorder'): # Test if a list or set has only zero elements iszero:=proc(v) local i; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: M. van Hoeij`; if type(v,{list,set}) then for i in v do if not procname(i) then return false fi od; true else i:=traperror(normal(v)); evalb(i=0 or not ( testeq(i)=false or traperror(simplify(evala(Normal(i))))<>0 ) ) fi end: #savelib('iszero'): # Compute root(A,n) ROOT:=proc(A,x,n) local B,i; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: M. van Hoeij`; B := traperror(sqrfree(A,x)); if B=lasterror then B:=factors(A) fi; root(B[1]*mul(i[1]^modp(i[2],n),i=B[2]),n,'symbolic') * mul(i[1]^floor(i[2]/n),i=B[2]) end: #savelib('ROOT'): #------------------------------------------------------------------------ is_remove_lower:=proc(A) option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved.`; iszero(A[1]) end: #savelib('is_remove_lower'): remove_lower := proc(A,x) local AA, i, j, n, y, k,L,Dx; option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved. \ Author: George Labahn and M. van Hoeij`; if not iszero(A[1]) then [] else n := nops(A); i := 1; k := 0; while i < n and iszero(A[i+1]) do i := i+1; od; if i = n-1 then AA := []; else # Addition, MVH April 2000: first attempt to integrate using # integrate_sols, because it is faster than int, and # usually produces nicer looking results. L:=add(A[i+1+j]*Dx^j,j=0..n-i-1); if not has(A,x) then k:=i # No need to integrate fi; while k lasterror and j::list then # integrate_sols succesful L:=j[1]; k:=k+1 else # integrate_sols unsuccesful, give up break fi od; AA := dsolve( coeff(L,Dx,0)*y(x) + add(coeff(L,Dx,j)*diff(y(x),x$j),j=1..n-i-1),y(x),'output'='basis'); fi; if i < n-1 then for j from 1 to i-k do # integrate_sols integrated k times, so the result needs # to be integrated another i-k times. AA := map(`dsolve/int`,AA,x); od; fi; [seq(x^(j-1),j=1..i), op(AA) ]; fi; end: #savelib('remove_lower'): #------------------------------------------------------------------------ # Pattern 1: L = symmetric_power(DF^2+a(x)*DF+b(x),m, [DF,x]) # where m=n-1. # # How to recognize pattern 1: # # L = DF^(m+1) + add(A.i * DF^i,i=0..m) where: # A.m = 1/2*a(x)*m*(m+1)*DF^m + # A.(m-1) = 1/24*m*(m-1)*(3*m+2)*(m+1)*a(x)^2+1/6*m*(m+1)*(m+2)*b(x) # +1/6*m*(m+1)*diff(a(x),x)*(m-1) # # The formula for A.m can easily be proved using the fact that # A.m=-normal(diff(Y,x)/Y) where # Y:=det(wronskian([seq(y1(x)^i*y2(x)^(m-i),i=0..m)],x)) # # The formula for A.(m-1) could be proved using the wronskian( # [y(x), seq(y1(x)^i*y2(x)^(m-i),i=0..m)],x), but it can be proved more # easily using Theorem 1 in the ISSAC'97 paper of Bronstein/Mulders/Weil. # # George Labahn, Mark van Hoeij, 1997. # Pattern 1. is_sympow_o2 := proc(A,x,v) local i, a, DF, b, m, compare, y; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Authors: George Labahn and M. van Hoeij`; m := nops(A) - 2; a := Normalizer( 2/m/(m+1)*A[m+1] ); b := Normalizer( 1/4*(-3*m^4*a^2-2*m^3*a^2-4*m^3*diff(a,x)+3*m^2*a^2 +2*m*a^2+4*m*diff(a,x)+24*A[m])/m/(m^2+3*m+2)); compare := symmetric_power(DF^2 + a*DF + b,m,[DF,x]); v:=[m, diff(y(x),x,x)+a*diff(y(x),x)+b*y(x)=0, y(x)]; iszero([seq(coeff(compare,DF,i)-A[i+1],i=0..m-2)]) end: #savelib('is_sympow_o2'): solve_sympow_o2:=proc(v) local i,ans; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: George Labahn and M. van Hoeij`; ans:=dsolve(v[2],v[3],'output'='basis'); if nops(ans) = 2 then [seq( ans[1]^i*ans[2]^(v[1]-i),i=0..v[1])] else [] fi end: #savelib('solve_sympow_o2'): #------------------------------------------------------------------------ # Pattern 2: L = LCLM(DF-RootOf(_Z^n+a(x)),groundfield=[]) # # How to recognize pattern 2: # # L = DF^n + add(A.i * DF^i,i=0..n-1) where: # A.0 = a(x) # A.(n-1)=-(n-1)/2*diff(a(x),x)/a(x) # If n>2 then A.(n-2) = # (n-1)*(n-2)*(7*n-1)/n/24 * (diff(a(x),x)/a(x))^2 # - (n-1)*(n-2)/6 * diff(a(x),x,x)/a(x); # # Notes: 1) This is a generalization of plus_min_expsols. # 2) If R = RootOf(_Z^n-a(x)) is a reducible RootOf then L is a # reducible equation. If DFactor has already been tried then # we know that R is an irreducible RootOf. # # George Labahn, Mark van Hoeij, June 2000. # Pattern 2. is_LCLM_nth_root := proc(A,x,v) local i,n,DF,compare,R; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Authors: G. Labahn and M. van Hoeij`; n:=nops(A)-1; if has(A[1],x) and iszero(A[n]+(n-1)/2*diff(A[1],x)/A[1]) and (n<3 or iszero( (n-1)*(n-2)*(7*n-1)/n/24 * (diff(A[1],x)/A[1])^2 - (n-1)*(n-2)/6 * diff(A[1],x,x)/A[1] - A[n-1])) and frontend(irreduc,[frontend(primpart,[DF^n+A[1],DF])]) then R:=RootOf(_Z^n+A[1]); if n>3 then # A[1], A[n] and A[n-1] have already been checked, so more # checks are only necessary when n>3 compare:=DEtools['LCLM'](DF-R,[DF,x],'groundfield'=[args]) fi; if (n<4 or iszero([seq(coeff(compare,DF,i)-A[i+1],i=0..n)])) then v:=[R,A[1],n]; return true fi fi; false end: #savelib('is_LCLM_nth_root'): solve_LCLM_nth_root:=proc(v,x) local ans,i,B; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Authors: G. Labahn and M. van Hoeij`; ans := EXP_INT(v[1],x); # = exp(int(R,x)) with exp(ln(..)) simplified. [seq(subs(v[1]=subs(B=v[2],convert(RootOf(_Z^v[3]+B,'index'=i) ,radical)),ans),i=1..v[3])] end: #savelib('solve_LCLM_nth_root'): #------------------------------------------------------------------------ # Pattern 3: L = symmetric_product( DF^2+a(x)-b(x),DF^2+a(x)+b(x) ) # # How to recognize pattern 3: # # L = DF^4 - B(x)*DF^3 + 4*a(x)*DF^2 # + (-4*B(x)*a(x)+6*diff(a(x),x))*DF # + 4*b(x)^2 - 2*B(x)*diff(a(x),x) + 2*diff(a(x),x,x) # # where B(x)=diff(b(x),x)/b(x) # so b(x)=exp(Int(B(x),x)) # # Additional properties: # Let L1 = symmetric_power(DF^2+a(x)-b(x),2) # L2 = symmetric_power(DF^2+a(x)+b(x),2) # and M1 = symmetric_power(L1,2), and same for M2. # Then # 1) symmetric_power(L,2) = symmetric_product(M1,M2) and # in particular the order of this is < 10. # # 2) exterior_power(L,2) = LCLM(M1, M2) # This statement is slightly different when the coeff of # DF^1 in DF^2+a(x)-b(x), DF^2+a(x)+b(x) is non-zero. # # Can also write L as: # # L = DF^4 + a3(x)*DF^3 + a2(x)*DF^2 # + (3/2*diff(a2(x),x)+a2(x)*a3(x))*DF # + DESol({-2*diff(a0(x),x)+diff(a2(x),x,x,x) # +diff(a3(x),x)*diff(a2(x),x)+3*a3(x)*diff(a2(x),x,x) # -4*a3(x)*a0(x)+2*a3(x)^2*diff(a2(x),x)},{a0(x)}) # # Mark van Hoeij, June 2000. # Pattern 3. is_sympr_o2_no_DF1 := proc(A,x,v) local y,bx2,bx,I1,I2; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: M. van Hoeij`; if nops(A)=5 and A[1]<>0 and iszero(3*diff(A[3],x)/2+A[3]*A[4]-A[2]) then bx2:=Normalizer((A[1]-diff(A[3],x,x)/2-2*A[4]*diff(A[3],x)/4)/4); # b(x)^2 if bx2<>0 and iszero(diff(bx2,x)/bx2 + 2*A[4]) then bx:=ROOT(bx2,x,2); I1,I2:=Normalizer(A[3]/4-bx), Normalizer(A[3]/4+bx); v:=[diff(y(x),x,x)+I1*y(x), diff(y(x),x,x)+I2*y(x), y(x)]; return true fi fi; false end: #savelib('is_sympr_o2_no_DF1'): solve_sympr_o2_no_DF1 := proc(v) local i,j,ans1,ans2; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: M. van Hoeij`; ans1 := dsolve(v[1],v[3],'output'='basis'); if not type(ans1,list) then ans1:=[ans1] fi; ans2 := dsolve(v[2],v[3],'output'='basis'); if not type(ans2,list) then ans2:=[ans2] fi; [seq(seq(i*j,i=ans1),j=ans2)] end: #savelib('solve_sympr_o2_no_DF1'): #------------------------------------------------------------------------ # Pattern 4: L = LCLM(DF^2+a(x)*DF+b(x)+RootOf(_Z^2+c(x)),groundfield=[]) # # How to recognize pattern 4: # # L = DF^4 + add(A[i+1] * DF^i,i=0..3) # and a few things must hold. First: exterior_power(L,2) has order 5 # instead of 6, which holds when the following is zero: # (4*A[3]-6*diff(A[4],x)-A[4]^2)*A[4]-8*A[2]+8*diff(A[3],x)-4*diff(A[4],x,x) # # If this is so, we first assume a(x)=0. Note that this implies that # coeff(exterior_power(L,2),DF,0) = 0 # From coeff(L,DF,3) we can determine C(x)=diff(c(x),x)/c(x) if we # assume a(x)=0, and b(x) can then be determined from coeff(L,DF,2). # We don't need to look at coeff(L,DF,1) because the condition # degree(exterior_power(L,2),DF)=5 leaves only one possible value for # coeff(L,DF,1) as a function of the other coeffs. Now from coeff(L,DF,0) # we can find c(x) in terms of b(x) and C(x). The check is now that # C(x) is the logarithmic derivative of c(x). If this holds, then the # assumptions that L is such LCLM and that a(x)=0 are correct. # If not, then we're going to try to find a(x), do a transformation # to make it zero and apply recursion. # Now coeff(L,DF,3) equals 2*a(x)-C(x). We can clear that term by # sending DF to DF-a(x)/2+C(x)/4 which moves # LCLM(DF^2+a(x)*DF+... # LCLM(DF^2+C(x)/2*DF+... # so we see that c(x)^(-1/2) = exp(int(-C(x)/2)) is a double solution # of the 2nd exterior power (1 of these causes the order to drop from # 6 to 5, but the other solution is still there, and it is the square # root of a rational function, so we're only looking for exponents # in 1/2 * integer). So we calculate a solution of the 2nd exterior # power that is a square root of a rational function, if it exists. # This gives c(x)^(-1/2), which lets us find C(x), and then a(x). # # Mark van Hoeij, June 2000. # Pattern 4. is_LCLM_o2_same_DF1 := proc(A,x,V) local dif,A4_1,A4_2,A4_3,A4_4,A3_1,A3_2,A3_3,A1_1,Cx,v,ax,bx,cx,Ls,R,I1,I2, DF, i; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: M. van Hoeij`; if nops(A)=5 and iszero( (4*A[3]-6*diff(A[4],x)-A[4]^2)*A[4] -8*A[2]+8*diff(A[3],x)-4*diff(A[4],x,x) ) then # Let c(x)=exp(Int(C(x),x)) Cx:=-A[4]; bx:=Normalizer(A[3]/2+diff(Cx,x)/4-1/8*Cx^2); # Now 2*diff(bx,x)-bx*Cx-A[2] equals zero. v:=diff(bx,x); cx:=Normalizer(A[1]+v*Cx-bx^2-diff(v,x)+bx/2*diff(Cx,x)-bx/4*Cx^2); if cx=0 then # cx can be zero, e.g. on mult( DF^2+x*DF+x^5, DF^2+x*DF+x^5 ) # but this reducible operator is not an LCLM. return false elif iszero(cx*Cx-diff(cx,x)) then # Now ax=0 R:=ROOT(-cx,x,2); I1,I2:=DF^2+Normalizer(bx+R), DF^2+Normalizer(bx-R); if nargs>3 and args[-2]=`2nd try` then I1,I2:=op(map(symmetric_product,[I1,I2],DF+args[-1],[DF,x])); fi; V:=[map(DEtools['diffop2de'],[I1,I2],y(x),[DF,x]), y(x)]; return true elif nargs=3 and type(A,list(ratpoly(anything,x))) then dif:=proc(a,x) Normalizer(diff(a,x)) end: A4_1:=dif(A[4],x); A4_2:=dif(A4_1,x); A4_3:=dif(A4_2,x); A4_4:=dif(A4_3,x); A3_1:=dif(A[3],x); A3_2:=dif(A3_1,x); A3_3:=dif(A3_2,x); A1_1:=dif(A[1],x); v:=map(Normalizer,[ -A4_4-A[3]*A4_2+A[3]*A3_1-2*A1_1+A3_3+(-A4_2-1/2*A3_1)*A4_1+(-A4_3+1/2* A3_2-1/2*A4_1*A[3]+(-1/4*A3_1+1/8*A[4]*A4_1)*A[4])*A[4], -4*A[1]+3*A3_2-3/4*A4_1^2-7/2*A4_3+(A3_1-11/4*A4_2+1/16*A[4]^3)*A[4]+ (-1/2*A[4]^2-2*A4_1+A[3])*A[3], 3*A3_1-9/4*A[4]*A4_1-9/2*A4_2, -3/4*A[4]^2-3*A4_1+2*A[3], 0,1]); v:=`DEtools/Expsols`(v,0,x,`use Int`,`no algext`,radical,denom,2); if v<>[] then v:=sort(v,(i,j) -> evalb(length(i) evalb(length(i)0 then userinfo(2,dsolve, `Multiplying solutions by`,exp(Int(Bx,x))); v:=symmetric_product(DF-Bx,add(A[j+1]*DF^j,j=0..4),[DF,x]); if is_sympr_o2_no_DF1([seq(coeff(v,DF,j),j=0..4)],x,'w') then V:=[op(w),-Bx]; return true fi fi od fi; false end: #savelib('is_sympr_o2'): solve_sympr_o2:=proc(v,x) local ans,i,j; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Author: M. van Hoeij`; ans:=solve_sympr_o2_no_DF1(v); j:=EXP_INT(v[4],x); [seq(j*i,i=ans)] end: #savelib('solve_sympr_o2'): # # ---> to_normal_form(A,x) # # If none of the routines in higherorder produce an answer and if # the second highest coefficient is nonzero then transform equation # to one with second highest coefficient zero and try some routines # again. # # Note: # # This trick is already used in dsolve in to_const_coeff where one # checks to see if ode transforms as above to constant coefficient # or Euler type equations. # # Note: # # 1) do not try symmetric power as previous call would have obtained answer # # 2) only try LCLM methods if transform changes ode from non-rational # function coefficients to one having rational functions coefficients. # # George Labahn, Mark Van Hoeij 2000 # to_normal_form := proc(A,x) local ord, ode, newA, ans, NF, DF, i, u, is_rational; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved. \ Authors: G. Labahn and Mark van Hoeij`; ord := nops(A)-1; ode := add(A[i+1]*DF^i,i=0..ord); ode := symmetric_product(ode, DF-A[ord]/ord, [DF,x]); newA:= [seq(coeff(ode,DF,i),i=0..ord)]; is_rational := evalb(type(A,list(ratpoly(anything,x)))); ans := higherorder( newA, x, `second try`, is_rational ); if ans <> [] then NF := EXP_INT(-A[ord]/ord,x); [seq(u*NF, u=ans)] else [] fi; end: #savelib('to_normal_form'):