share/contrib/maxima-odesolve/kovacicODE.mac

   1 /* kovacic.mac
   2
   3   Solve second order linear ODEs with Liouvillian solutions using Kovacic' algorithm
   4
   5   References:
   6   [1] Carolyn J. Smith A discussion and implementation of Kovacic' algorithm,
   7       MSC thesis university of waterloo, 1984
   8       https://cs.uwaterloo.ca/research/tr/1984/CS-84-35.pdf
   9   [2] B.D. Saunders An implementation of Kovacic's algorithm for solving second order linear homogeneous differential equations
  10   SYMSAC '81 Proceedings of the fourth ACM symposium on Symbolic and algebraic computation
  11   Pages 105-108
  12   [3] Jerald J. Kovacic, An algorithm for solving second order linear homogeneous differential equations, Journal of Symbolic Computation, v.2 n.1, p.3-43, March 1986
  13
  14
  15   Copyright (C) 2014 Nijso Beishuizen
  16
  17   This program is free software; you can redistribute it and/or modify
  18   it under the terms of the GNU General Public License as published by
  19   the Free Software Foundation; either version 2 of the License, or
  20   (at your option) any later version.
  21
  22   This program is distributed in the hope that it will be useful,
  23   but WITHOUT ANY WARRANTY; without even the implied warranty of
  24   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  25   GNU General Public License for more details.
  26
  27   You should have received a copy of the GNU General Public License
  28   along with this program; if not, write to the Free Software
  29   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  30 */
  31
  32 /*****************************************************************************************************/
  33 put('kovacic,001,'version)$
  34 /*****************************************************************************************************/
  35
  36 /* ----- print everything with flag lower than DEBUGFLAG ----- */
  37 DEBUGFLAG:1$
  38
  39 /* load absimp to get rid of abs(a) when a>0 */
  40 /*load(absimp)$*/
  41
  42 /* NOTES */
  43 /* remember: one of the necessary conditions might be too strict, see kovacic, compare with Saunders !!!*/
  44 /* we should check for linear independence by checking that the wronskian has nonzero determinant*/
  45 /*     y1 y2  */
  46 /* W = y1' y2'*/
  47 /*            */
  48 /* (det(W))  = (y1y2' -y1'y2')  */
  49 /* (det(W))' = (y1y2' -y1'y2')' */
  50 /* so detW is in C */
  51 /* 1a. check the paper of Ulmer & Weil */
  52 /* 1b. check paper of Man, combining with Prelle-Singer algorithm */
  53 /* 2: check if linearizable by lie point transformation */
  54 /* 3. check if linearizable by generalized Sundman transformation */
  55 /* 4. check if linearizable by euler-Liouville transformation */
  56 /* 5. check if linearizable by cartan equivalence      */
  57 /* 6. check if transformable from linear second order without liouvillian  */
  58 /*      _solutions to linear second order with liouvillian _solutions  */
  59 /*      for this we need to have a mapping to a set of standard linear equations */
  60 /*      for instance, how does y"=xy transform into y"=0 ?*/
  61
  62
  63 /*****************************************************************************************************/
  64 /* ----- Check if y(x)=_solution is a _solution of (differential) equation                     ----- */
  65 /*****************************************************************************************************/
  66 isSolution(_solution,_equation):=block([_eq,_subs1],
  67   /* assumes equation contains "=" */
  68   _eq:rhs(_equation)-lhs(_equation),
  69   /* assumes _solution has the form y(x)=expression */
  70   _subs1:subst(_solution,_equation),
  71   return(ratsimp(ev(_subs1,diff)))
  72 )$
  73 /*****************************************************************************************************/
  74
  75
  76
  77 /*****************************************************************************************************/
  78 /* ----- Check if expr is an equation (containing a "=" )                                      ----- */
  79 /*****************************************************************************************************/
  80 isEquation(_expr) := block([],
  81   if freeof("=",_expr) then return(false) else return(true)
  82 )$
  83
  84 /* rules for linear second order equation: */
  85 /*   y'' = b(x)*y'+c(x)*y + d(x)                 */
  86
  87 /*****************************************************************************************************/
  88 matchdeclare ([_b,_c,_d],freeof(_dy))$
  89 matchdeclare ([_c1,_c2],freeof(_x))$
  90
  91 /* match for linear second order equation - can match dependencies:
  92   depends(y,x)
  93   ode2_linear(a*diff(y,x),y,x)
  94 as well as regular _variables:
  95   ode2_linear(a*diff(y(x),x),y(x),x)    */
  96 defmatch (ode2_linear, _b*'diff(_dy,_x)+_c*_dy+_d,_dy,_x)$
  97 defmatch (linear, _c1*_x+_c2,_x)$
  98
  99
 100
 101 /*****************************************************************************************************/
 102 /*****************************************************************************************************/
 103 kovacicODE(_expr,_y1,_x1) := block([_y,_x,_var,_ode_order,_ddy,_ode,_phi,_b,_c,_d,_s,_t,_r,_ord_inf,_l,_m,_n,_j,_trial,_squo,_srem,_tcont,_nu],
 104 /*****************************************************************************************************/
 105 /*****************************************************************************************************/
 106   /* ----- substitute our own local _variables ----- */
 107   _expr:subst(_y,_y1,_expr),
 108   _expr:subst(_x,_x1,_expr),
 109
 110   /* what do we know about _x */
 111   /*print(facts(_x1)), */
 112
 113   _var:_x,
 114
 115   dprint(0,"equation: ",_expr,", degree = ",derivdegree(_expr,_y,_x)),
 116   /* ----- first check if it is an equation (we only check for equal sign) ----- */
 117   if isEquation(_expr)=false then (
 118     dprint(0,_expr, " is not an equation."),
 119     return(false)
 120   ),
 121
 122   /* ----- find the order of the ODE ----- */
 123   _ode_order : derivdegree(_expr,_y,_x),
 124   if (_ode_order=2) then (
 125     dprint(0,"Second order ODE found.")
 126   )
 127   else (
 128     dprint(0,_expr," is not a second order ODE."),
 129     return(false)
 130   ),
 131
 132
 133   /* ----- make sure the equation is in the form y'' = F(x,y,y') ----- */
 134   _ddy: 'diff(_y,_x,2),
 135   /* we only take the first _solution, in general, loop over the _solutions! */
 136   _ode: solve(_expr,_ddy)[1],
 137   _phi : rhs(_ode),
 138   if (lhs(_ode)=_ddy) and (freeof(_ddy,_phi)) then (
 139     dprint(1,"ODE: ",string(_ode))
 140   )
 141   else (
 142     dprint(0,"could not separate second order differential operator ",_ode),
 143     return(false)
 144   ),
 145
 146
 147   /*print("is ode2 linear: ",ode2_linear(_phi,_y,_x)),*/
 148   /* ----- check if the equation is linear ----- */
 149   if (ode2_linear(_phi,_y,_x)=false) then (
 150     dprint(0,"ODE is not linear!"),
 151     return(false)
 152   )
 153   else (
 154     dprint(1,"ODE is linear: ",_b,_c,_d),
 155     _b:-_b,
 156     _c:-_c,
 157     _d:-_d
 158   ),
 159
 160   /*print(ev(_phi,_y=0)),*/
 161
 162   /* ----- check if equation is homogeneous ----- */
 163   if (ev(_phi,_y = 0)=0) then (
 164     dprint(1,"ODE is homogeneous:")
 165   )
 166   else (
 167     dprint(1,"equation is nonhomogeneous, first trying homogeneous case")
 168     /*_d : 0*/
 169     /* could be transformed to third order homogeneous ode l(y)=0 */
 170     /* L(y)=b  --->   l(y) = b*(L(y))' -b'*L(y) = 0  */
 171     /*return(false)*/
 172   ),
 173
 174   /* ----- Transform to normal form y''= (s/t)*y, ----- */
 175   _s : 2*diff(_b,_x) + _b*_b - 4*_c,
 176   _t : 4,
 177   dprint(5,"s = ",_s),
 178   dprint(5,"t = ",_t),
 179
 180
 181   /* remember, maple's normal(in expanded form) is like maxima's ratsimp */
 182   _r : ratsimp(_s/_t),
 183   dprint(5,"normal = ",_r),
 184   _s : ratexpand(num(_r)),
 185   _t : ratexpand(denom(_r)),
 186   dprint(5,"s=",_s),
 187   dprint(5,"t=",_t),
 188
 189   /* ----- determine the quotient and remainder of _s/_t for the main _variable _x ----- */
 190   /* catch this because divide is not exactly polynomial long division and the quotient is not uniquely determined */
 191   if freeof(_x,_s) and freeof(_x,_t) and (_t#0) then (
 192     [_squo,_srem] : [_s/_t,0]
 193   )
 194   else (
 195     [_squo,_srem] : divide(_s,_t,_x),
 196     _squo:ratexpand(_squo), /* we need to do a ratexpand because we want to get the coefficient correctly*/
 197     _srem:ratexpand(_srem) /* we need to do a ratexpand because we want to get the coefficient correctly*/
 198   ),
 199   dprint(5,"quo = ",_squo),
 200   dprint(5,"rem = ",_srem),
 201
 202   /* ----- compute a square free factorization ----- */
 203   [_sdec,_tcont] : sqfree(_t,_var),
 204   dprint(5,"_sdec square free = ",_sdec),
 205
 206   /* this is the order*/
 207   _m : length(_sdec),
 208   dprint(5,"m = ",_m),
 209
 210   _t : _tcont,
 211
 212
 213   /* TODO: can we replace part(_sdec,i) with _sdec[i] ? */
 214   for _i:1 thru _m do (
 215     _t : _t*part(_sdec,_i)^_i
 216   ),
 217
 218   if _m>0 then (
 219     _t1 : part(_sdec,1)
 220   )
 221   else (
 222     _t1 : 1
 223   ),
 224
 225   if _m>1 then (
 226     _t2: part(_sdec,2)
 227   )
 228   else (
 229    _t2 : 1
 230   ),
 231
 232   dprint(5,"_t2 = ",_t2),
 233   dprint(5,"t = ",_t),
 234   dprint(5,"_t1 = ",_t1),
 235   dprint(5,"_t2 = ",_t2),
 236   dprint(5,"s = ",_s),
 237
 238   /* this does not give the same results as kovacic examples in the 1986 paper: is Smith wrong? */
 239   _ord_inf : hipow(ratexpand(_t),_var) - hipow(ratexpand(_s),_var),
 240   dprint(1,"order at infinity = ",_ord_inf),
 241
 242   _listl : [],
 243
 244   _oddti : true,
 245   for _i from 3 step 2 thru _m do (
 246     if (part(_sdec,_i)#1) then (
 247       dprint(2,"_oddti is false, it cannot be case 1"),
 248       _oddti : false
 249     )
 250   ),
 251
 252   if _oddti and (featurep(_ord_inf/2,integer) or (_ord_inf>2)) then (
 253     dprint(2,"case 1 applies"),
 254     _listl : append(_listl,[1])
 255   ),
 256
 257   if not _oddti or (_t2 # 1) then (
 258     dprint(2,"case 2 applies"),
 259     _listl : append(_listl,[2])
 260   ),
 261
 262   if (_m <= 2) and (_ord_inf >= 2) then (
 263     _listl : append(_listl,[4,6,12])
 264   ),
 265
 266   dprint(1,"list of cases that apply: L=",_listl),
 267
 268   if length(_listl)=0 then (
 269     dprint(0,"No Liouvillian solutions exist!"),
 270     return(false)
 271   ),
 272
 273
 274 /* ************************************************* */
 275   dprint(0, "end of preliminaries..."),
 276 /* ************************************************* */
 277
 278   /* STEP 1 PART (a) */
 279   _dfix : (min(_ord_inf,2) - hipow(ratexpand(_t),_var) - 3*hipow(ratexpand(_t1),_var))/4,
 280   _thetafix : ratsimp((diff(_t,_var)/_t + 3*diff(_t1,_var)/_t1)/4),
 281   dprint(5,"_dfix=",_dfix),
 282   dprint(5,"_thetafix=",_thetafix),
 283
 284   /* STEP 1 PART (b) */
 285   /* Poles of order 2: find the roots c1,...,c_k2 of _t2 */
 286   /* first, get a list of all the unique roots */
 287   _rlist2 : rootz(_t2,_var),
 288   dprint(5,"_rlist2 = ",_rlist2),
 289
 290   /* find the leading coefficient of _t2 */
 291   /*lcoeff : ratcoef(ratexpand(_t2),_var,hipow(ratexpand(_t2),_var)),*/
 292   _t2 : _t2/lcoeff(_t2,_var),
 293   dprint(5,"_t2=",_t2),
 294
 295   /* k2 is the number of roots of _t2*/
 296   _k2:length(_rlist2),
 297   dprint(5,"_k2 = ",_k2),
 298
 299   /* loop over all roots */
 300   for _i:1 thru _k2 do (
 301     dprint(5,"i=",_i),
 302
 303     trest : _t/_t2^2 * product((_var-_rlist2[_j])^2,_j,1,_i-1) * product((_var-_rlist2[_j])^2,_j,_i+1,_k2),
 304     dprint(5,"srem  = ",_srem),
 305     dprint(5,"trest = ",trest),
 306     dprint(5,"var = ",_var),
 307     dprint(5,"_rlist2 = ",_rlist2[_i]),
 308
 309     _sol:undetcoeff(_srem,trest,_var,_rlist2[_i],2,2),
 310     dprint(5,"_sol:",_sol),
 311     _dd[_i] : ratsimp(1+4*_sol)^(1/2),
 312     _theta[_i] : ratsimp(_dd[_i]/(_var-_rlist2[_i])),
 313     dprint(5,"d1 = ",_dd[_i]),
 314     dprint(5,"_theta1 = ",_theta[_i])
 315   ),
 316
 317   dprint(5,"m = ",_m),
 318
 319   /* high order poles */
 320   /* STEP 3 PART (c) */
 321   _k1 : _k2,
 322
 323 if member(1,_listl) then (
 324   dprint(1,"1 is in L so we calculate higher order poles"),
 325   for _i:4 thru _m step 2 do (
 326     dprint(5,"_i = ",_i, "/",_m),
 327     L_i : _sdec[_i],
 328     dprint(5,"L_i=",L_i),
 329     _rlisthigher : rootz(L_i,_var),
 330     /*lcoeff : ratcoef(ratexpand(L_i),_var,hipow(ratexpand(L_i),_var)),*/
 331     dprint(5,"sdec = ",_sdec),
 332     dprint(5,"_i-1:",_i-1),
 333     dprint(5,"_m:",_m),
 334     dprint(5,"part1 = ",create_list(_sdec[ii],ii,1,_i-1)),
 335     dprint(5,"part2 = ",L_i/lcoeff(L_i,_var)),
 336     dprint(5,"part3 = ",create_list(_sdec[ii],ii,_i+1,_m)),
 337     /* */
 338     _sdec : append(create_list(_sdec[ii],ii,1,_i-1), [L_i/lcoeff(L_i,_var)], create_list(_sdec[ii],ii,_i+1,_m)),
 339     dprint(5,"sdec = ",_sdec),
 340     _nu : _i/2,
 341     /* we could do for rt in _rlisthigher ??? */
 342     for _j from 1 thru length(_rlisthigher) do (
 343       _k1 : _k1 + 1,
 344       rt : _rlisthigher[_j], /* !!! */
 345       dprint(5,_t,_sdec[_i],_i),
 346       trest : _t/_sdec[_i]^_i * product((_var-_rlisthigher[_l])^_i,_l,1,_j-1) * product((_var-_rlisthigher[_l])^_i,_l,_j+1,length(_rlisthigher)),
 347       dprint(5,_srem,trest),
 348       /* method of undetermined coefficients */
 349       _sol:undetcoeff(_srem,trest,_var,rt,2*_nu,2*_nu),
 350       if (_sol=0) then (
 351         dprint(0,"No Liouvillian solutions exist"),
 352         return(false)
 353       ),
 354       dprint(0, "_sol = ",_sol),
 355       _ac[_nu] : ratsimp(_sol^(1/2)),
 356       for _k from _nu-1 thru 2 step -1 do (
 357         _ac[_k] : _vtemp,
 358         result : sum(_ac[_l]*_ac[_nu+_k-_l],_l,_k,_nu),
 359
 360         /* can we do _ac[k]=undetcoeff directly because result is already solved?*/
 361         undet : undetcoeff(_srem,trest,_var,rt,2*_nu,_k+_nu),
 362         dprint(5,"result:",result),
 363         dprint(5,"undet:",undet),
 364         _ac[_k] : rhs(solve(result=undet,_vtemp)[1]),
 365         dprint(5,"_ac[_k]=",_ac[_k])
 366       ),
 367
 368       /* sum of _ac[k]*_ac[_nu+1-k] */
 369       result : sum(_ac[_k]*_ac[_nu+1-_k],_k,2,_nu-1),
 370
 371       _dd[_k1] : (undetcoeff(_srem,trest,_var,rt,2*_nu,_nu+1)-result)/_ac[_nu],
 372
 373       result : sum(_ac[_k]/(_var-rt)^_k,_k,2,_nu),
 374
 375       _theta[_k1] : 2*result + _dd[_k1]/(_var-rt)
 376
 377     )
 378   )
 379
 380 ),
 381 dprint(3,"_dd = ",_dd),
 382
 383 if _ord_inf>2 then (
 384   dprint(1,"_ord_inf > 2"),
 385   _dd[0]:1,
 386   _theta[0]:0
 387 ) else
 388 if _ord_inf = 2 then (
 389   dprint(1,"_ord_inf=2"),
 390   dprint(5,"s=",_s),
 391   dprint(5,"t=",_t),
 392
 393   _sol:lcoeff(_s,_var)/lcoeff(_t,_var),
 394   _dd[0]:ratsimp((1+4*_sol)^(1/2)),
 395   _theta[0]:0
 396 )
 397 else if member(1,_listl) then (
 398   dprint(1,"1 is in L"),
 399   _nu:(-_ord_inf)/2,
 400   dprint(5,"_nu = ",_nu),
 401   dprint(5,"_squo = ",_squo),
 402   dprint(5,"_var = ",_var),
 403   dprint(5,"coeff=",coeff(_squo,_var,2*_nu)),          /* should we just use lcoeff? it uses ratexpand..*/
 404   _ac[_nu]:ratsimp(coeff(_squo,_var,2*_nu)^(1/2)),
 405   dprint(5,"_ac[_nu]=",_ac[_nu]),
 406   for _i from _nu-1 thru 0 step -1 do (
 407     _ac[_i]:_vtemp,
 408     dprint(5,"_i=",_i),
 409     dprint(5,"_ac[_i]=",_ac[_i]),
 410
 411     _temp_var : sum(_ac[_j]*_ac[_i+_nu-_j],_j,_i,_nu),
 412
 413     dprint(5,"_temp_var = ",_temp_var),
 414     dprint(5,"coeff = ",coeff(_squo,_var,_i+_nu)),
 415     _ac[_i]:rhs(solve(_temp_var=coeff(_squo,_var,_i+_nu),_vtemp)[1]),
 416     dprint(5,"_ac[_i]=",_ac[_i])
 417   ),
 418   /* note: ac is not a list, we have just defined ac[1],ac[2]*/
 419   dprint(5,"1.ac = ",_ac),
 420
 421   dprint(5,"_nu = ",_nu),
 422   _temp_var : sum(_ac[_l]*_ac[_nu-_l-1],_l,0,_nu-1),
 423   dprint(5,"tempvar = ",_temp_var),
 424
 425   if _nu=0 then (
 426     _t : ratexpand(_t),
 427     _srem : ratexpand(_srem),
 428     _squo : ratexpand(_squo),
 429     dprint(5,"lcoeff=",lcoeff(_t,_var)),
 430     dprint(5,"degree = ",hipow(_t,_var)),
 431     dprint(5,"coeff=",coeff(_srem,_var,hipow(_t,_var)-1)),
 432     _aa: coeff(_srem,_var,hipow(_t,_var)-1),
 433     _bb: lcoeff(_t,_var),
 434     if (_bb#0) then (
 435       _tmp1:_aa/_bb - _temp_var
 436     )
 437     else (
 438       /* or should we simply bail out when the denominator=0?*/
 439       _tmp1: -_temp_var
 440     )
 441   )
 442   else (
 443     _tmp1:coeff(_squo,_var,_nu-1) - _temp_var
 444   ),
 445   dprint(5,"_tmp1=",_tmp1),
 446
 447   /* the denominator can be zero - we catch it by saying that d[0] should be zero then */
 448   /* but is this true?*/
 449   if (_tmp1=0 or _ac[_nu]=0) then (
 450     _dd[0]:0
 451   )
 452   else (
 453     _dd[0]:_tmp1/_ac[_nu]
 454   ),
 455   dprint(0, "_dd[0] = ",_dd[0]),
 456
 457   _tmp2 : sum(_ac[_l]*_var^_l,_l,0,_nu),
 458
 459   _theta[0]:2*_tmp2
 460 )
 461 else (
 462   dprint(2,"otherwise"),
 463   _dd[0]:0,
 464   _theta[0]:0
 465 ),
 466
 467
 468 dprint(3,"d0 = ",_dd[0]),
 469 dprint(3,"_theta0 = ",_theta[0]),
 470
 471 /* ************************************************* */
 472 /* *****  step 2 - form trial d's and _theta's ***** */
 473 /* ************************************************* */
 474 dprint(1,"Now entering step 2: forming trial ds and _thetas, length= ",length(_listl)),
 475
 476 _solfound:0,
 477
 478 for i from 1 thru length(_listl) do (
 479   dprint(1,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"),
 480   dprint(1,"!!!   CASE ",_listl[i], " ",i,"/",length(_listl),"      !!!"),
 481   dprint(1,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"),
 482
 483   _n:_listl[i],
 484
 485   /* Saunders */
 486   if _n = 1 then (
 487     _m : _k1
 488   )
 489   else (
 490     _m : _k2
 491   ),
 492
 493   dprint(3,"m = ",_m),
 494
 495   if (_n=2) and (_ord_inf<2) then (
 496     _dd[0]:0,
 497     _theta[0]:0
 498   ),
 499
 500   /* s = (-n/2,-n/2+1,...n/2)*/
 501   for _j from 0 thru _m do (
 502     dprint(5,"entering loop!"),
 503     /* -n/2, -n/2, -n/2, -n/2, ... */
 504     sq[_j]:-1/2 * _n
 505   ),
 506   /* in case of n=1 (case 1) this is (-1/2, +1/2) */
 507   _solution:[],
 508   alls:false,
 509   /* for all sequences s={-n/2, -n/2+1, ..., n/2}*/
 510
 511   /* preliminary: d0 should be */
 512   _dslist:[],
 513   while (alls # true) do (
 514     dprint(5,"entering do loop"),
 515
 516     /*tezt : sq[0]*d[0] - sum(sq[_l]*d[_l],_l,1,_m),*/
 517     tezt : sq[0]*_dd[0],
 518     for _l from 1 thru _m do (
 519       tezt : tezt - sq[_l]*_dd[_l]
 520     ),
 521
 522     ds:ratsimp(_n*_dfix+tezt),
 523     if not member(ds,_dslist) then (
 524       _dslist:append(_dslist,[ds]),
 525       dprint(3,"dslist:",_dslist),
 526       dprint(3,"ds = ",ds, ", should be integer and >=0",featurep(ds,integer), " ",is(ds>=0)),
 527       if (integerp(ds)=true) then (
 528         if (ds>=0) then (
 529         dprint(1,"ds is a positive or zero integer"),
 530         tmp : sum(sq[_l]*_theta[_l],_l,0,_m),
 531         dprint(5,"tmp : ",tmp),
 532         _thetas:ratsimp(_n*_thetafix+tmp),
 533         dprint(5,"_thetas = ",_thetas),
 534         /* step3 - determine polynomial P if possible and hence omega and _solution */
 535         dprint(5,"step 3"),
 536         _soln:step3(_n,ds,_thetas,_s/_t,_var),
 537         dprint(5,"finished step 3"),
 538         _soln : radcan(_soln), /* to simplify */
 539         dprint(3,"_soln = ",_soln,length([_soln])),
 540
 541         if _soln#false then (
 542           dprint(0, "finding sol2"),
 543           ratio:_b,
 544
 545           /*gamma_expand:true, */ /* to get rid of incomplete gamma in saunders ex. 1 */
 546           /* actually, better to work with gamma to prevent complex logarithms like kamke196 */
 547           /* kamke 2.111: we should get exp(-x)/x here */
 548
 549           _soln1:exp(integrate(-1/2*ratio,_var))*_soln, /* list [soln] necessary?*/
 550           dprint(5,"_soln1 = ",_soln1),
 551
 552           /* kill the impaginary part so we don't have to cancel out all the imaginary bits */
 553           _soln1:trigsimp(radcan(_soln1)),
 554           if not(freeof(%i,_soln1)) then (
 555             _soln1 : trigsimp(realpart(_soln1))
 556           ),
 557           _soln1:factor(_soln1),
 558           _soln1:ev(_soln1,integrate),
 559
 560           /* logarc gets rid of arcsinh stuff inside exponentials */
 561           /* sometimes this increases the complexity of the solution, so do a check afterwards */
 562           _soln11:radcan(logarc(_soln1)),
 563
 564           /* if logarc has lead to a better solution, then use it */
 565           _nrOps11 : nrOps(_soln11),
 566           _nrOps1 : nrOps(_soln1),
 567
 568           /* use the radcan-logarc solution when the expression contains less operations
 569              or less types of operations */
 570           if ((length(_nrOps11) <= length(_nrOps1)) or (length(unique(_nrOps11))<=length(unique(_nrOps1))) ) then (
 571             _soln1 : _soln11
 572           ),
 573
 574           _soln1:radcan(_soln1),
 575
 576           dprint(3,"x. _soln1 = ",_soln1),
 577
 578           /*soln2 :radcan(soln1*integrate((exp(-integrate(ratio,_var))/(soln1*soln1)),_var)),*/
 579           _soln2:-integrate(ratio,_var),
 580           _soln2:exp(_soln2)/(_soln1*_soln1),
 581           _soln22:radcan(logarc(_soln2)),
 582
 583           /* if logarc has lead to a better solution, then use it */
 584           _nrOps22 : nrOps(_soln22),
 585           _nrOps2 : nrOps(_soln2),
 586
 587           /* use the radcan-logarc solution when the expression contains less operations
 588              or less types of operations */
 589           if ((length(_nrOps22) <= length(_nrOps2)) or (length(unique(_nrOps22))<=length(unique(_nrOps2))) ) then (
 590             _soln2 : _soln22
 591           ),
 592
 593           dprint(3,"x1. _soln2 = ",_soln2,", ",_var),
 594
 595           _soln2:integrate(_soln2,_var),
 596           dprint(3,"x2. _soln2 = ",_soln2),
 597           _soln2:trigsimp(radcan(_soln1*_soln2)),
 598           dprint(3,"x3. _soln2 = ",_soln2),
 599           if not(freeof(%i,_soln2)) then (
 600           dprint(3,"x3a. _soln2 = ",_soln2),
 601
 602             _soln2:trigsimp(realpart(_soln2))  /* take the real part and drop the imaginary part */
 603           ),
 604           dprint(3,"x4. _soln2 = ",_soln2),
 605           _soln2 : radcan(_soln2), /*  to further simplify - we still have exp(2^(3/2)*x-sqrt(2)*x) here */
 606           dprint(3,"x5. _soln2 = ",_soln2),
 607           _soln2 : ev(_soln2,nouns),
 608           dprint(3,"length=",length([_soln])),
 609           if length([_soln]) = 1 then (
 610             dprint(3,"_soln 1"),
 611             _solution: [ratsimp(_soln1),ratsimp(_soln2)], /* realpart because we can have imaginary solutions for y'' + y = 0  and we don't want to deal with the trouble of cancelling them out elegantly by itself, so we just delete them*/
 612             dprint(3,"homogeneous solution:",_solution),
 613             _solfound : 1,
 614             return(true)
 615           )
 616           else (
 617             dprint(3,"_soln 2"),
 618             _solution: [ratsimp(_soln1),ratsimp(_soln2)],op(2,[_soln]),
 619             _solfound:1,
 620             return(true)
 621           )
 622         )
 623         else (
 624           dprint(0,"No Liouvillian solutions exist"),
 625           _solution: false
 626         ),
 627
 628       dprint(5,"alls:",alls)
 629       )
 630     )
 631   ), /* end of dsloop*/
 632
 633   _jj : _m-1,
 634   dprint(5,"sq_jj = ",sq),
 635   for _j from _m thru 0 step -1 do (
 636       /*print("j = ",_j,"/",_m), */
 637       if (sq[_j] = (1/2 * _n)) then (
 638         _jj : _jj - 1,
 639         sq[_j]: -1/2 * _n
 640       )
 641       else (
 642         _jj : 0,
 643         sq[_j] : sq[_j] + 1,
 644         /* exit (the for loop) */
 645         return(true)
 646       )
 647     ),
 648
 649     if (_jj < 0) then (
 650       dprint(0,"No Liouvillian solutions exist"),
 651       _solution: false,
 652       alls:true
 653     )
 654
 655   ),
 656
 657   if _solfound=1 then (
 658     return(true)
 659   )
 660 ),
 661
 662
 663   if (alls=false) then (
 664   /* ----- nonhomogeneous _solution -----*/
 665   dprint(1,"checking if a nonhomogeneous part exists"),
 666   _P:0,
 667   if (_d#0) then (
 668     /* 2 independent _solutions of the homogeneous equation */
 669     _s : _solution,
 670     dprint(1,"nonhomogeneous part found - looking for particular solution"),
 671     /* -------------------------------------------------- */
 672     /* ----- METHOD 1                                     */
 673     /* -------------------------------------------------- */
 674     _E : exp(integrate(_b,_var)),
 675     _PP : integrate(_E*_s[1]*_d,_var),
 676     /* we need radcan because of the quotient by zero bug */
 677     _PP : (_PP/(_E*_s[1]*_s[1])),
 678     _PP : radcan(_PP),
 679     /* there is an error somewhere - see kamke 184*/
 680     _PP : -_s[1]*integrate(_PP,_var),
 681
 682     /* perform the simplification after we have compared the two solutions */
 683     /* we do this because simplification can lead to expression swell, like for */
 684     /* example kamke 2.234 */
 685     /* so we hope that simplification of the smallest term will still lead to */
 686     /* a simpler expression than simplification of the largest term */
 687     /* -------------------------------------------------- */
 688     /* ----- METHOD 2                                     */
 689     /* -------------------------------------------------- */
 690     /* Determinant of the Wronskian */
 691     _W : _s[2]*diff(_s[1],_var) - diff(_s[2],_var)*_s[1],
 692     _W : ratexpand(_W),
 693     _W : radcan(_W),
 694     _W : trigsimp(_W),
 695     _W : trigreduce(_W),
 696     _W : ratsimp(_W),
 697     dprint(5,"5. wronskian : ", _W),
 698     dprint(1,"nonhomogeneous part : ", _d),
 699     /* particular _solution of the ODE */
 700     _P : _s[2]*integrate(_s[1]*_d/_W,_var) - _s[1]*integrate(_s[2]*_d/_W,_var),
 701
 702     dprint(1,"0. particular solution : ", _P),
 703     /*_P : ratexpand(_P),*/
 704     /*_P : trigexpand(_P),*/
 705     /*
 706     _P : trigreduce(_P),
 707     dprint(0, "1.  particular solution : ", _P),
 708     _P : trigsimp(_P),
 709     dprint(0, "2. particular solution : ", _P),
 710     _P : radcan(_P),
 711     dprint(0, "3. particular solution : ", _P),
 712     */
 713     /* substitute back the original _variables */
 714     /* -------------------------------------------------- */
 715
 716
 717     /*_P : subst([_x=_x1],_P),*/
 718     /* use method one, because it supposedly can give less complicated particular solutions */
 719     /*print("x=",_x),*/
 720     /*print("x1=",_x1),*/
 721
 722     /*
 723     dprint(0, "complexity of method 1 for particular solution: ",length(nrOps(_PP))),
 724     dprint(0, "complexity of method 2 for particular solution: ",length(nrOps(_P))),
 725     */
 726     /* if method 1 has less operation than method 2, use method 1 */
 727     if (length(nrOps(_PP)) <= length(nrOps(_P)) ) then (
 728       _P : subst([_x=_x1],_PP)
 729     )
 730     else (
 731       _P : subst([_x=_x1],_P)
 732     ),
 733
 734     _P : trigreduce(_P),
 735     _P : trigsimp(_P),
 736     _P : radcan(_P),
 737     dprint(5,"3. particular solution : ", _P)
 738   ),
 739
 740
 741   dprint(1,"substituting original variables back"),
 742   _solution[1]:subst([_x=_x1,_y=_y1], _solution[1]),
 743   _solution[2]:subst([_x=_x1,_y=_y1], _solution[2]),
 744
 745 /*
 746   _solution[1]:trigsimp(trigreduce(_solution[1])),
 747   _solution[2]:trigsimp(trigreduce(_solution[2])),
 748   _solution[1]:ev(_solution[1],nouns),
 749   _solution[2]:ev(_solution[2],nouns),
 750   _solution[1]:trigsimp(trigreduce(_solution[1])),
 751   _solution[2]:trigsimp(trigreduce(_solution[2])),
 752 */
 753
 754   _C1 : constant_factors(factor(_solution[1])),
 755   _C2 : constant_factors(factor(_solution[2])),
 756   _solution[1] : ratsimp(_solution[1]/_C1), /* ratsimp necessary for kamke 2.286*/
 757   _solution[2] : ratsimp(_solution[2]/_C2),
 758   _solution: [_y1 = _solution[1]*%k1 + _solution[2]*%k2 + _P],
 759
 760 /*_solution:[subst([_x=_x1,_y=_y1], _y = _solution[1]*%k1 + _solution[2]*%k2 + _P)],*/
 761
 762
 763
 764
 765   dprint(1,"11d:",string(_solution))
 766   ),
 767
 768   dprint(1,"solution:",string(_solution)),
 769   /* also return the method used to solve the ode */
 770   method : 'kovacic,
 771   return(_solution)
 772 )$
 773
 774
 775
 776
 777 /*****************************************************************************************************/
 778 step3(_n,_d,_theta,rhs1,_var):=block([_p,_listv,_i,_a,_a0,_a1,_a2,_a3,_a4,_a5,_pr,_sete,_soln,_trial,_w],
 779 /*****************************************************************************************************/
 780 /* ************************************************* */
 781 /* *****   step 3 -                            ***** */
 782 /* ************************************************* */
 783
 784   /*# form P in terms of undetermined coefficients a_i */
 785   /*# P = a-{d}*x^{d} + a_{d-1}*x^{d-1} + ... + a_0*/
 786   dprint(3,"_var=",_var),
 787   dprint(3,"degree of polynomial = ",_d),
 788   dprint(3,"n=",_n),
 789   dprint(3,"_theta=",_theta),
 790   dprint(3,"rhs1=",rhs1),
 791
 792   /* start with highest order*/
 793   _p:_var^_d,
 794   dprint(3,"P = ",_p),
 795
 796   /* [a_d-1, a_d-2, ..., a_2, a_1, a_0] */
 797   _listv : makelist(concat(_a,_d-1 -_i),_i,_d-1,0,-1),
 798
 799   _p : _p + sum(_listv[_i]*_var^(_i-1),_i,1,length(_listv)),
 800
 801   dprint(3,"_listv=",_listv),
 802
 803   /*# generate recursive relations P_i */
 804   _pr[_n]: -_p,
 805   for _i from _n thru 0 step -1 do (
 806     _pr[_i-1]: ratsimp(-diff(_pr[_i],_var) - _theta * _pr[_i] - (_n-_i)*(_i+1)*rhs1*_pr[_i+1])
 807   ),
 808
 809   dprint(3,"P-1 = ",_pr[-1]),
 810   _trial:ratexpand(num(ratsimp(_pr[-1]))),
 811   dprint(3,"_trial = ",_trial),
 812
 813   if _trial # 0 then (
 814     _sete : [], /* */
 815     dprint(3,"low degree = ",lopow(_trial,_var)),
 816     dprint(3,"high degree = ",hipow(_trial,_var)),
 817
 818     for _i from lopow(_trial,_var) thru hipow(_trial,_var) do (
 819       lset:coeff(_trial,_var,_i),
 820       dprint(3,"lset:",lset),
 821       if (lset # []) then (
 822         _sete:endcons(lset,_sete) /* create list as [sete,lset] */
 823       )
 824     ),
 825     dprint(5,"sete = ",_sete),
 826     dprint(5,"_listv = ",_listv),
 827
 828     if(length(_listv)=0) then (
 829       dprint(0,"No Liouvillian solutions exist"),
 830       return(false)
 831     ),
 832
 833     dprint(5,"a0 = ",_a0),
 834     _soln:solve(_sete,_listv),
 835
 836     dprint(1,"_soln=",_soln),
 837     /*if (_soln[1]=[] or length(_soln)=0 ) then (*/
 838     if (length(_soln)=0 ) then (
 839       dprint(1,"no solutions"),
 840       return(false)
 841     ),
 842
 843     /* for more than 1 element, the solution has two brackets: [[1,2,3,...]]*/
 844     if (length(_listv)#1) then (
 845       _soln:_soln[1]
 846     ),
 847     /*_soln:_soln[1],*/
 848     dprint(5,"determining ais"),
 849     /* ai = _soln[d-_i] */
 850     /* finally, we are getting some stuff here for kamke 2.129 */
 851     for _i from _d-1 thru 0 step -1 do(
 852       map(":",[concat(_a,_i)],[rhs(_soln[_d-_i])])
 853     )
 854   ),
 855   dprint(5,"a0 = ",_a0),
 856   dprint(5,"P = ",_p),
 857   _trial : sum(_pr[_i]*_w^_i/(_n-_i)!,_i,0,_n),
 858
 859
 860   dprint(3,"_trial = ",string(_trial)),
 861   dprint(3,"w = ",_w),
 862
 863   _solution:solve(ev(_trial),_w),
 864
 865   if _solution=[] then (
 866     dprint(2,"w is empty: ",_w),
 867     dprint(2,"NOTE: we still need to correct the return value!"),
 868     return(exp(int(_w,_var)),_trial=0)
 869   ),
 870   dprint(1,"we have a genuine _solution:",_solution),
 871
 872   _w:ratsimp(factor(_solution[1])), /* we need factor, it's like bringing it in normal form */
 873
 874 /* this call to ratexpand causes radcan to behave strangely (doesn't simplify certain things) later on */
 875 /*  dprint(0, "expanded = ",ratexpand(_w)),*/
 876   dprint(3,"w=",_w),
 877   sol:exp(integrate(rhs(_w),_var)),
 878   dprint(3,"sol=",sol),
 879   return(sol)
 880 )$
 881
 882
 883
 884 /*****************************************************************************************************/
 885 /*****************************************************************************************************/
 886 lcoeff(_t2,_var):=block([],
 887  return(ratcoef(ratexpand(_t2),_var,hipow(ratexpand(_t2),_var)))
 888 )$
 889
 890
 891 /*****************************************************************************************************/
 892 /* ----- calculate the roots and expand the multiple roots                                     ----- */
 893 /*****************************************************************************************************/
 894 rootz(_expr,_var):=block([_sol,_mult,_rlist2,_j,_i],
 895   /* note: this has the potential to fail */
 896   _sol:solve(_expr,_var),
 897   /* get the multiplicities of the roots*/
 898   _mult:multiplicities,
 899   /* create a list of roots including multiple roots */
 900   _rlist2: flatten( makelist( makelist(rhs(_sol[_j]),_i,1,_mult[_j]),_j,1,length(_sol))),
 901   return(_rlist2)
 902 )$
 903
 904
 905 /*****************************************************************************************************/
 906 /* ----- determine coefficient of a factor in a partial fraction expansion                     ----- */
 907 /*****************************************************************************************************/
 908 undetcoeff(_num,_rden,_var,_root,_m,_ex):=block([_k,_p,_er],
 909   dprint(0, "undetcoeff:",_m-_ex,_rden,_num),
 910   _k : _m - _ex,
 911   _p : _num/_rden,
 912
 913   /* differentiate p with respect to _var k times */
 914   _p : ratsimp(diff(_p,_var,_k)),
 915   dprint(5,"_p:",_p),
 916   /* catch errors, it now returns empty when an error occurs (e.g. when we have non-polynomial expressions)  */
 917   _er: errcatch(subst(_root,_var,_p)/_k!),
 918   if (_er=[]) then (
 919     _er:0
 920   ) else (
 921     _er:_er[1]
 922   ),
 923   dprint(0, "_er:",_er),
 924   return(_er)
 925 )$
 926
 927
 928 /*****************************************************************************************************/
 929 /* ----- Transform second order ode y"=f(x)y'+g(x)y to normal form y"=F(x)y ----- */
 930 /*****************************************************************************************************/
 931 Normalform(_phi,_y,_x) := block([_f,_g,_df,_F],
 932   _f : coeff(ratexpand(_phi),_y),
 933   _g : coeff(ratexpand(_phi),'diff(_y,_x)),
 934   _df : diff(_f,_x),
 935   _F: radcan(ratsimp(_g - _f*_f/4 - _df/2)),
 936   return(_F)
 937 )$
 938
 939
 940 /*****************************************************************************************************/
 941 /* ----- transforms _expr from normal form solution to actual solution using rhs of the ode    ----- */
 942 /*****************************************************************************************************/
 943 NormalSolutionToActualSolution(_expr,_phi) :=block([_f,_g,_df,_sol],
 944   _f : coeff(ratexpand(_phi),_y),
 945   _g : coeff(ratexpand(_phi),'diff(_y,_x)),
 946   _df : diff(_f,_x),
 947   _sol : radcan(ratsimp(_expr*exp(integrate(_f,_x)/2))),
 948   return(_sol)
 949 )$
 950
 951
 952 /*****************************************************************************************************/
 953 /* ----- return the primary part of the polynomial                                              -----*/
 954 /*****************************************************************************************************/
 955 primpart(_expr,_x):=block([],
 956   return(_expr/content(_expr,_x))
 957 )$
 958
 959
 960 /*****************************************************************************************************/
 961 /* ----- Squar-free decomposition (Yun, On square free decomposition algorithms)               ----- */
 962 /*****************************************************************************************************/
 963 sqfree(_expr,_var):=block([_i,_signp,_tc,_tlist1,_c,_d,_n,_cc,_re,_pp,_cont,_g,_PPP],
 964
 965   /* if polynomial is a constant, then return the constant */
 966   if freeof(_var,_expr) then (
 967     dprint(0, "polynomial is constant"),
 968     return[_expr,0]
 969   ),
 970
 971   _expr:ratexpand(_expr),
 972
 973   /* get the highest power of the expression, sorted by leading terms first, then get the operator */
 974   _n:hipow(_expr,_var),
 975
 976   if atom(_expr) then (
 977     _signp:1
 978   )
 979   else (
 980     _cc:coeff(_expr,_var,_n),
 981     if atom(_cc) then (
 982       _signp:1
 983     )
 984     else (
 985       if op(_cc)="-" then (
 986         _signp:-1
 987       )
 988       else (
 989         _signp:1
 990       )
 991     )
 992   ),
 993
 994   _re: ratexpand(_expr)/_signp,
 995   /* --- */
 996   _pp: primpart(_re,_var)[1],
 997   _tc:content(_re,_var)[1],
 998   _cont:_tc*_signp,
 999   /*--- */
1000
1001   _tlist1:[],
1002
1003   _re:ratexpand(diff(_pp,_var)),
1004   _g:gcd(_pp,_re),
1005   [_PPP,_c,_d]:[_g,ratsimp(_pp/_g),ratsimp(_re/_g)],
1006
1007   while _c#1 do (
1008     _re: ratexpand(_d-diff(_c,_var)),
1009     _g:gcd(_c,_re),
1010     [_PPP,_c,_d]:[_g,ratsimp(_c/_g),ratsimp(_re/_g)],
1011     _tlist1 : append(_tlist1,[_PPP])
1012   ),
1013
1014   return([_tlist1,_cont])
1015 )$
1016
1017
1018 /*****************************************************************************************************/
1019 /* ----- simple method of finding the constant factor in front of an equation                  ----- */
1020 /*****************************************************************************************************/
1021 constant_factors(_expr) := block([inflag:true],
1022            if not mapatom(_expr) and op(_expr)="*"
1023            then xreduce("*",listify(subset(setify(args(_expr)),'constantp)))
1024            else 1
1025 )$
1026
1027
1028
1029 /*****************************************************************************************************/
1030 /* ----- calculates the number of operators in the expression ----- */
1031 /* ----- this is a simple measure of complexity               ----- */
1032 /*****************************************************************************************************/
1033 nrOps(expression):=
1034  block( [ ],
1035         dprint(0, "expression:",expression),
1036         allOpsPriv (expression, [])
1037        )$
1038
1039 allOpsPriv(expression, opList) :=
1040  block ( [x, args, newList],
1041         /* if expression is an atom, then we return opList */
1042         if atom(expression)
1043            then opList
1044            else (
1045               x:    op(expression),
1046               args: args(expression),
1047               /* add the operators to the opList */
1048               newList: cons(x, opList),
1049               for arg in args do
1050                 /* also expand all the subexpressions and count the operators */
1051                 newList: allOpsPriv(arg, newList),
1052               newList
1053              )
1054         )$
1055
1056
1057 /*****************************************************************************************************/
1058 /* ----- dprint(flag,expr) is debuggingprint. It only prints expr when flag<DEBUGFLAG          ----- */
1059 /*****************************************************************************************************/
1060 dprint(flag,[expr])::= if flag < DEBUGFLAG then buildq ([expr], print (splice (expr)));