axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / RDETRS.spad

)abbrev package RDETRS TranscendentalRischDESystem
++ Author: Manuel Bronstein
++ Date Created: 17 August 1992
++ Date Last Updated: 3 February 1994
++ References
++ Bron90 The Transcendental Risch Differential Equation
++ Description:
++ Risch differential equation system, transcendental case.

TranscendentalRischDESystem(F, UP) : SIG == CODE where
  F  : Join(Field, CharacteristicZero, RetractableTo Integer)
  UP : UnivariatePolynomialCategory F
 
  N   ==> NonNegativeInteger
  Z   ==> Integer
  RF  ==> Fraction UP
  V   ==> Vector UP
  U   ==> Union(List UP, "failed")
  REC ==> Record(z1:UP, z2:UP, r1:UP, r2:UP)
 
  SIG ==> with

    monomRDEsys : (RF, RF, RF, UP -> UP) -> _
                   Union(Record(a:UP, b:RF, h:UP, c1:RF, c2:RF, t:UP),"failed")
      ++ monomRDEsys(f,g1,g2,D) returns \spad{[A, B, H, C1, C2, T]} such that
      ++ \spad{(y1', y2') + ((0, -f), (f, 0)) (y1,y2) = (g1,g2)} has a solution
      ++ if and only if \spad{y1 = Q1 / T, y2 = Q2 / T},
      ++ where \spad{B,C1,C2,Q1,Q2} have no normal poles and satisfy
      ++ A \spad{(Q1', Q2') + ((H, -B), (B, H)) (Q1,Q2) = (C1,C2)}
      ++ D is the derivation to use.

    baseRDEsys : (RF, RF, RF)   -> Union(List RF, "failed")
      ++ baseRDEsys(f, g1, g2) returns fractions \spad{y_1.y_2} such that
      ++ \spad{(y1', y2') + ((0, -f), (f, 0)) (y1,y2) = (g1,g2)}
      ++ if \spad{y_1,y_2} exist, "failed" otherwise.
 
  CODE ==> add

    import MonomialExtensionTools(F, UP)
    import SmithNormalForm(UP, V, V, Matrix UP)
 
    diophant: (UP, UP, UP, UP, UP) -> Union(REC, "failed")
    getBound: (UP, UP, UP, UP, UP) -> Z
    SPDEsys : (UP, UP, UP, UP, UP, Z, UP -> UP, (F, F, F, UP, UP, Z) -> U) -> U
    DSPDEsys: (F, UP, UP, UP, UP, Z, UP -> UP) -> U
    DSPDEmix: (UP, UP, F, F, N, Z, F) -> U
    DSPDEhdom: (UP, UP, F, F, N, Z) -> U
    DSPDEbdom: (UP, UP, F, F, N, Z) -> U
    DSPDEsys0: (F, UP, UP, UP, UP, F, F, Z, UP -> UP, (UP,UP,F,F,N) -> U) -> U
 
-- reduces (y1', y2') + ((0, -f), (f, 0)) (y1,y2) = (g1,g2) to
-- A (Q1', Q2') + ((H, -B), (B, H)) (Q1,Q2) = (C1,C2), Q1 = y1 T, Q2 = y2 T
-- where A and H are polynomials, and B,C1,C2,Q1 and Q2 have no normal poles.
-- assumes that f is weakly normalized (no finite cancellation)
    monomRDEsys(f, g1, g2, derivation) ==
      gg := gcd(d := normalDenom(f, derivation),
               e := lcm(normalDenom(g1,derivation),normalDenom(g2,derivation)))
      tt := (gcd(e, differentiate e) exquo gcd(gg,differentiate gg))::UP
      (u := ((tt * (aa := d * tt)) exquo e)) case "failed" => "failed"
      [aa, tt * d * f, - d * derivation tt, u::UP * e * g1, u::UP * e * g2, tt]
 
-- solve (y1', y2') + ((0, -f), (f, 0)) (y1,y2) = (g1,g2) for y1,y2 in RF
-- assumes that f is weakly normalized (no finite cancellation) and nonzero
-- base case: F' = 0
    baseRDEsys(f, g1, g2) ==
      zero? f => error "baseRDEsys: f must be nonzero"
      zero? g1 and zero? g2 => [0, 0]
      (u := monomRDEsys(f, g1, g2, differentiate)) case "failed" => "failed"
      n := getBound(u.a, bb := retract(u.b), u.h,
                    cc1 := retract(u.c1), cc2 := retract(u.c2))
      (v := SPDEsys(u.a, bb, u.h, cc1, cc2, n, differentiate,
                   (z1,z2,z3,z4,z5,z6) +->
                    DSPDEsys(z1, z2::UP, z3::UP, z4, z5, z6, differentiate)))
                          case "failed" => "failed"
      l := v::List(UP)
      [first(l) / u.t, second(l) / u.t]
 
-- solve
--   D1 = A Z1 + B R1 - C R2
--   D2 = A Z2 + C R1 + B R2
-- (D1,D2) = ((A, 0, B, -C), (0, A, C, B)) (Z1, Z2, R1, R2)
-- for R1, R2 with degree(Ri) < degree(A)
-- assumes (A,B,C) = (1) and A and C are nonzero
    diophant(a, b, c, d1, d2) ==
      (u := diophantineSystem(matrix [[a,0,b,-c], [0,a,c,b]],
                          vector [d1,d2]).particular) case "failed" => "failed"
      v := u::V
      qr1 := divide(v 3, a)
      qr2 := divide(v 4, a)
      [v.1 + b * qr1.quotient - c * qr2.quotient,
       v.2 + c * qr1.quotient + b * qr2.quotient, qr1.remainder, qr2.remainder]
 
-- solve
-- A (Q1', Q2') + ((H, -B), (B, H)) (Q1,Q2) = (C1,C2)
-- for polynomials Q1 and Q2 with degree <= n
-- A and B are nonzero
-- cancellation at infinity is possible
    SPDEsys(a, b, h, c1, c2, n, derivation, degradation) ==
      zero? c1 and zero? c2 => [0, 0]
      n < 0 => "failed"
      g := gcd(a, gcd(b, h))
      ((u1 := c1 exquo g) case "failed") or
        ((u2 := c2 exquo g) case "failed") => "failed"
      a := (a exquo g)::UP
      b := (b exquo g)::UP
      h := (h exquo g)::UP
      c1 := u1::UP
      c2 := u2::UP
      (da := degree a) > 0 =>
        (u := diophant(a, h, b, c1, c2)) case "failed" => "failed"
        rec := u::REC
        v := SPDEsys(a, b, h + derivation a, rec.z1 - derivation(rec.r1),
                    rec.z2 - derivation(rec.r2),n-da::Z,derivation,degradation)
        v case "failed" => "failed"
        l := v::List(UP)
        [a * first(l) + rec.r1, a * second(l) + rec.r2]
      ra := retract(a)@F
      ((rb := retractIfCan(b)@Union(F, "failed")) case "failed") or
        ((rh := retractIfCan(h)@Union(F, "failed")) case "failed") =>
                                DSPDEsys(ra, b, h, c1, c2, n, derivation)
      degradation(ra, rb::F, rh::F, c1, c2, n)
 
-- solve
-- a (Q1', Q2') + ((H, -B), (B, H)) (Q1,Q2) = (C1,C2)
-- for polynomials Q1 and Q2 with degree <= n
-- a and B are nonzero, either B or H has positive degree
-- cancellation at infinity is not possible
    DSPDEsys(a, b, h, c1, c2, n, derivation) ==
      bb := degree(b)::Z
      hh:Z :=
        zero? h => 0
        degree(h)::Z
      lb := leadingCoefficient b
      lh := leadingCoefficient h
      bb < hh =>
        DSPDEsys0(a,b,h,c1,c2,lb,lh,n,derivation,
                  (z1,z2,z3,z4,z5) +-> DSPDEhdom(z1,z2,z3,z4,z5,hh))
      bb > hh =>
        DSPDEsys0(a,b,h,c1,c2,lb,lh,n,derivation,
                  (z1,z2,z3,z4,z5) +-> DSPDEbdom(z1,z2,z3,z4,z5,bb))
      det := lb * lb + lh * lh
      DSPDEsys0(a,b,h,c1,c2,lb,lh,n,derivation,
                (z1,z2,z3,z4,z5) +-> DSPDEmix(z1,z2,z3,z4,z5,bb,det))
 
    DSPDEsys0(a, b, h, c1, c2, lb, lh, n, derivation, getlc) ==
      ans1 := ans2 := 0::UP
      repeat
        zero? c1 and zero? c2 => return [ans1, ans2]
        n < 0 or (u:= getlc(c1,c2,lb,lh,n::N)) case "failed" => return "failed"
        lq := u::List(UP)
        q1 := first lq
        q2 := second lq
        c1 := c1 - a * derivation(q1) - h * q1 + b * q2
        c2 := c2 - a * derivation(q2) - b * q1 - h * q2
        n := n - 1
        ans1 := ans1 + q1
        ans2 := ans2 + q2
 
    DSPDEmix(c1, c2, lb, lh, n, d, det) ==
      rh1:F :=
        zero? c1 => 0
        (d1 := degree(c1)::Z - d) < n => 0
        d1 > n => return "failed"
        leadingCoefficient c1
      rh2:F :=
        zero? c2 => 0
        (d2 := degree(c2)::Z - d) < n => 0
        d2 > n => return "failed"
        leadingCoefficient c2
      q1 := (rh1 * lh + rh2 * lb) / det
      q2 := (rh2 * lh - rh1 * lb) / det
      [monomial(q1, n), monomial(q2, n)]
 
 
    DSPDEhdom(c1, c2, lb, lh, n, d) ==
      q1:UP :=
        zero? c1 => 0
        (d1 := degree(c1)::Z - d) < n => 0
        d1 > n => return "failed"
        monomial(leadingCoefficient(c1) / lh, n)
      q2:UP :=
        zero? c2 => 0
        (d2 := degree(c2)::Z - d) < n => 0
        d2 > n => return "failed"
        monomial(leadingCoefficient(c2) / lh, n)
      [q1, q2]
 
    DSPDEbdom(c1, c2, lb, lh, n, d) ==
      q1:UP :=
        zero? c2 => 0
        (d2 := degree(c2)::Z - d) < n => 0
        d2 > n => return "failed"
        monomial(leadingCoefficient(c2) / lb, n)
      q2:UP :=
        zero? c1 => 0
        (d1 := degree(c1)::Z - d) < n => 0
        d1 > n => return "failed"
        monomial(- leadingCoefficient(c1) / lb, n)
      [q1, q2]
 
-- return a common bound on the degrees of a solution of
-- A (Q1', Q2') + ((H, -B), (B, H)) (Q1,Q2) = (C1,C2), Q1 = y1 T, Q2 = y2 T
-- cancellation at infinity is possible
-- a and b are nonzero
-- base case: F' = 0
    getBound(a, b, h, c1, c2) ==
      da := (degree a)::Z
      dc :=
        zero? c1 => degree(c2)::Z
        zero? c2 => degree(c1)::Z
        max(degree c1, degree c2)::Z
      hh:Z :=
        zero? h => 0
        degree(h)::Z
      db := max(hh, bb := degree(b)::Z)
      da < db + 1 => dc - db
      da > db + 1 => max(0, dc - da + 1)
      bb >= hh => dc - db
      (n := retractIfCan(leadingCoefficient(h) / leadingCoefficient(a)
                      )@Union(Z, "failed")) case Z => max(n::Z, dc - db)
      dc - db