share/translators/m2mj/pmint.mpl

   1 # The Poor Man's Integrator, a parallel integration heuristic
   2 # Version 1.1 ---  May 10, 2005  (c) M.Bronstein and INRIA 2004-2005
   3 pmint := proc(f,x)
   4         local ff, si, li, lin, lout, ld, q, d, l, vars, dx, ls, fint, lc;
   5         ff := eval(convert(f, tan));            # convert trigs to tan
   6         si := select(proc(d) diff(d,x) <> 0 end, indets(ff));
   7         si := select(proc(d) diff(d,x) <> 0 end, indets(map(diff, si, x))) union si;
   8         li := [op(si)];         # list of terms in integrand and its derivative
   9         lin := [seq(d=`tools/genglobal`(x), d=li)];     # substitution terms->indets
  10         lout := [seq(rhs(d)=lhs(d), d=lin)];            # substitution indets->terms
  11         ld := subs(lin, map(diff, li, x));      # derivatives of all the terms
  12         q := lcm(seq(denom(d), d=ld));          # denominator of the total derivation
  13         l := [seq(normal(q * d), d=ld)];        # q * derivatives of all the terms
  14         vars := map(lhs, lout);
  15         dx := totalDerivation(vars, l);         # vector field dx = q * d/dx
  16         ls := [seq(getSpecial(d, lin), d=li)];  # list of known Darboux for dx
  17         fint := subs(lout, pmIntegrate(subs(lin, ff), dx, q, vars, ls));
  18         lc := select(proc(d) convert(d,string)[1]="_"; end, indets(fint, name));
  19         subs({seq(d = 0, d=lc minus si)}, fint);
  20 end;
  21
  22 getSpecial := proc(f, l) local p;               # return known Darboux polys
  23         p := op(0,f);
  24         if p = `tan` then [1+subs(l,f)^2, false];
  25         elif p = `tanh` then [1 + subs(l,f), false], [1 - subs(l,f), false];
  26         elif p = `LambertW` then [subs(l,f), true];
  27         else NULL; fi;
  28 end;
  29
  30 totalDerivation := proc(lv, ld)
  31         proc(f) local fp, i;
  32                 fp := 0; for i to nops(ld) do fp := fp + ld[i] * diff(f, lv[i]); od;
  33                 fp;
  34         end;
  35 end;
  36
  37 pmIntegrate := proc(f, d, q, vars)
  38         local ls, splq, s, ff, df, spl, cden, dg, monomials, cand, lunk, sol, i;
  39         if nargs = 5 then ls := args[5]; else ls := []; fi;
  40         splq := splitFactor(q, d);
  41         s := splq[1]; for i to nops(ls) do if ls[i][2] then s := s*ls[i][1]; fi; od;
  42         ff := normal(f); df := denom(ff); spl := splitFactor(df, d);
  43         cden := s * spl[1] * deflation(spl[2], d);
  44         dg := 1 + degree(s) +  max(degree(numer(ff)), degree(denom(ff)));
  45         monomials := [op(enumerateMonoms(vars, dg))];
  46         cand := add('_A'[i] * monomials[i], i = 1..nops(monomials)) / cden;
  47         lunk := { seq('_A'[i], i = 1..nops(monomials)) };
  48         sol:= tryIntegral(f, d, q, vars, cand, lunk, spl[1], spl[2], splq[1], ls, 0);
  49         if sol[1] then sol := tryIntegral(f, d, q, vars, cand, lunk,
  50                                                 spl[1], spl[2], splq[1], ls, I); fi;
  51         if sol[1] then Int(f); else sol[2]; fi;
  52 end;
  53
  54 tryIntegral := proc(f, d, q, vars, cand, lunk, l1, l2, l3, ls, K)
  55         local candlog, p, candidate, i, sol;
  56         candlog := [op({ myfactors(l1, K), myfactors(l2, K), myfactors(l3, K) }
  57                         union { seq(p[1], p=ls) })];
  58         candidate := cand + add('_B'[i] * log(candlog[i]), i = 1..nops(candlog));
  59         sol := solve({coeffs(numer(normal(f - d(candidate)/q)), {op(vars)})},
  60                                 lunk union { seq('_B'[i], i = 1..nops(candlog)) });
  61         [evalb(sol=NULL), subs(sol,candidate)];
  62 end;
  63
  64 myfactors := proc(p, K) local l, fact;
  65         if K = 0 then l := factors(p); else l := factors(p, K); fi;
  66         seq(fact[1], fact=l[2]);
  67 end;
  68
  69 enumerateMonoms := proc(vars, d) local n, x, i, v, s, w;
  70         n := nops(vars);
  71         if n = 0 then {1}; else
  72                 x := vars[n];
  73                 v := [seq(vars[i], i = 1..n-1)];
  74                 s := enumerateMonoms(v, d);
  75                 for i to d do s := s union {seq(x^i*w,w=enumerateMonoms(v,d-i))}; od;
  76                 s;
  77         fi;
  78 end;
  79
  80 splitFactor := proc(p, d) local si, x, c, q, spl, s, splh;
  81         si := select(proc(z) d(z) <> 0 end, indets(p,name));
  82         if si = {} then RETURN([1,p]) fi;
  83         x := si[1];
  84         c := content(p, x, 'q');
  85         spl := splitFactor(c, d);
  86         s := normal(gcd(q, d(q)) / gcd(q, diff(q, x)));
  87         if degree(s) = 0 then RETURN([spl[1], q * spl[2]]); fi;
  88         splh := splitFactor(normal(q / s), d);
  89         [spl[1] * splh[1] * s, spl[2] * splh[2]];
  90 end;
  91
  92 deflation := proc(p, d) local si, x, c, q;
  93         si := select(proc(z) d(z) <> 0 end, indets(p,name));
  94         if si = {} then RETURN(p) fi;
  95         x := si[1];
  96         c := content(p, x, 'q');
  97         deflation(c, d) * gcd(q, diff(q, x));
  98 end;