share/diffequations/ode2.mac

   1  /* -*- Mode: MACSYMA; Package: MAXIMA -*- */
   2
   3 /*  (c) Copyright 1981 Massachusetts Institute of Technology  */
   4
   5 /*  The Ordinary Differential Equation Solver.
   6     This package consists primarily of a set of routines taken from Moses'
   7     thesis and Boyce & DiPrima for solving O.D.E.s of 1st and 2nd order.
   8     The top-level routines are ODE2, IC1, IC2, and BC2.  */
   9
  10 eval_when(translate,
  11 declare_translated(boundtest,noteqn,nlxy,nly,nlx,xcc2,bessel2,euler2,
  12    pttest,exact2,cc2,genhom,solvebernoulli,solvehom,integfactor,exact,
  13    separable,solvelnr,solve1,linear2,reduce,hom2,pr2,varp,desimp,failure,
  14    ode1a,ftest,ode2a));
  15
  16 block(local(eq,yold,x),
  17 ode2(eq,yold,x):=block([derivsubst:false,ynew],
  18    subst(yold,ynew,ode2a(subst(ynew,yold,eq),ynew,x))))$
  19
  20 block(local(eq,y,x),
  21 ode2a(eq,y,x):=block([de,a1,a2,a3,a4,%q%],
  22    intfactor: false, method: 'none,
  23    if freeof('diff(y,x,2),eq)
  24      then if ftest(ode1a(eq,y,x)) then return(%q%) else return(false),
  25    if derivdegree(de: desimp(lhs(eq)-rhs(eq)),y,x) # 2
  26      then return(failure(msg1,eq)),
  27    a1: coeff(de,'diff(y,x,2)),
  28    a2: coeff(de,'diff(y,x)),
  29    a3: coeff(de,y),
  30    a4: expand(de - a1*'diff(y,x,2) - a2*'diff(y,x) - a3*y),
  31    if pr2(a1) and pr2(a2) and pr2(a3) and pr2(a4) and
  32       ftest(hom2(a1,a2,a3))
  33      then if a4=0 then return(%q%) else return(varp(%q%,-a4/a1)),
  34    if ftest(reduce(de)) then return(%q%) else return(false)))$
  35
  36 block(local(eq,y,x),
  37 ode1a(eq,y,x):=block([de,des], /* f, g, %q% */
  38    if derivdegree(de: expand(lhs(eq)-rhs(eq)),y,x) # 1
  39      then return(failure(msg1,eq)),
  40    if linear2(de,'diff(y,x)) = false then return(failure(msg2,eq)),
  41    des: desimp(de),
  42    de: solve1(des,'diff(y,x)),
  43    if ftest(solvelnr(de)) then return(%q%),
  44    if ftest(separable(de)) then return(%q%),
  45    if ftest(integfactor(%g%,%f%)) then return(exact(%q%*%g%,%q%*%f%)),
  46                                /* linear2 binds %f% and %g% */
  47    if linear2(des,'diff(y,x)) = false then return(failure(msg2,eq)),
  48    if ftest(integfactor(%g%,%f%)) then return(exact(%q%*%g%,%q%*%f%)),
  49    if ftest(solvehom(de)) then return(%q%),
  50    if ftest(solvebernoulli(de)) then return(%q%),
  51    if ftest(genhom(de)) then return(%q%) else return(false)))$
  52
  53 block(local(eq),
  54 desimp(eq):=block([inflag:true],
  55    eq: factor(eq),
  56    if atom(eq) or not(inpart(eq,0)="*") then return(expand(eq)),
  57    eq: map(lambda([u], if freeof(nounify('diff),u) then 1 else u), eq),
  58    return(expand(eq))))$
  59
  60 block(local(%f%),
  61 pr2(%f%):=freeof(y,'diff(y,x),'diff(y,x,2),%f%))$
  62
  63 block(local(call),
  64 ftest(call):=is((%q%: call) # false))$
  65
  66 block(local(eq,y),
  67 solve1(eq,y):=block([programmode:true],first(solve(eq,y))))$
  68
  69 /* linear2 tests for the form fx+%g% */
  70
  71 block(local(expr,x),
  72 linear2(expr,x):=block([],
  73    %f%: ratcoef(expr,x),
  74    if not(freeof(x,%f%)) then return(false),
  75    %g%: ratsimp(expr - %f%*x),
  76    return(freeof(x,%g%))))$
  77
  78 /*  variables used to denote constants: %C, %K1, %K2.
  79     METHOD denotes the method of solution.
  80     INTFACTOR denotes the integrating factor.
  81     ODEINDEX denotes the index for Bernoulli's method or for the genhom method.
  82     YP denotes the particular solution for the variation of parameters technique.  */
  83 \f
  84 /*  B&DiP, pp. 13-14  */
  85
  86 block(local(eq),
  87 solvelnr(eq):=block([%f%,%g%,w,%c],
  88    if linear2(rhs(eq),y) = false then return(false),
  89    w: %e^(integrate(%f%,x)),
  90    method: 'linear,
  91    return(y=w*(integrate(%g%/w,x)+%c))))$
  92
  93 /*  B&DiP, pp. 29-34  */
  94
  95 block(local(eq),
  96 separable(eq):=block([xpart:[],ypart:[],flag:false,inflag:true,%c],
  97    eq: factor(rhs(eq)),
  98    if atom(eq) or not(inpart(eq,0)="*") then eq: [eq],
  99    for u in eq do
 100       if freeof(x,u) then ypart: cons(u,ypart) else
 101       if freeof(y,u) then xpart: cons(u,xpart) else return(flag: true),
 102    if flag = true then return(false),
 103    if xpart = [] then xpart: 1 else xpart: apply("*",xpart),
 104    if ypart = [] then ypart: 1 else ypart: apply("*",ypart),
 105    method: 'separable,
 106    return(ratsimp(integrate(1/ypart,y)) = ratsimp(integrate(xpart,x)) + %c)))$
 107
 108 /*  B&DiP, pp. 34-41  */
 109
 110 block(local(m,n),
 111 integfactor(m,n):=block([b1,b2,dmdx,dmdy,dndx,dndy,dd,%e_to_numlog:true],
 112    dmdy: ratsimp(diff(m,y)),  dndx: ratsimp(diff(n,x)),
 113    if (dd: dmdy-dndx) = 0 then return(1),
 114    dmdx: ratsimp(diff(m,x)),  dndy: ratsimp(diff(n,y)),
 115    if dmdx-dndy=0 and dmdy+dndx=0 then return(1/(m^2 + n^2)),
 116    if freeof(y, (b1: ratsimp(dd/n))) then return(%e^(integrate(b1,x))),
 117    if freeof(x, (b2: ratsimp(dd/m)))
 118      then return(%e^(integrate(-b2,y))) else return(false)))$
 119
 120 block(local(m,n),
 121 exact(m,n):=block([a,ynew,%c],
 122    intfactor: subst(yold,ynew,%q%),
 123    a: integrate(ratsimp(m),x),
 124    method: 'exact,
 125    return(ratsimp(a + integrate(ratsimp(n-diff(a,y)),y)) = %c)))$
 126
 127 /*  B&DiP, pp. 43-44  */
 128
 129 block(local(eq),
 130 solvehom(eq):=block([qq,a1,a2,%c],
 131    a1: ratsimp(subst(x*qq,y,rhs(eq))),
 132    if not(freeof(x,a1)) then return(false),
 133    a2: ratsimp(subst(y/x,qq,integrate(1/(a1-qq),qq))),
 134    method: 'homogeneous,
 135    return(%c*x = %e^a2)))$
 136
 137 /*  B&DiP, p. 21, problem 15  */
 138
 139 block(local(eq),
 140 solvebernoulli(eq):=block([a1,a2,n,%c],
 141    a1: coeff(eq: expand(rhs(eq)),y,1),
 142    if not(freeof(y,a1)) then return(false),
 143    n: hipow(ratsimp(eq-a1*y),y),
 144    a2: coeff(eq,y,n),
 145    if not(freeof(y,a2)) or not(freeof(x,y,n)) or n=0
 146                         or not(eq = expand(a1*y + a2*y^n))
 147      then return(false),
 148    a1: integrate(a1,x),
 149    method: 'bernoulli, odeindex: n,
 150    return(y = %e^a1 * ((1-n)*integrate(a2*%e^((n-1)*a1),x) + %c) ^ (1/(1-n)))))$
 151 \f
 152 /*  Generalized homogeneous equation:  y' = y/x * H(yx^n)
 153         Reference:  Moses' thesis.  */
 154
 155 block(local(eq),
 156 genhom(eq):=block([%g%,u,n,a1,a2,a3,%c],
 157    %g%: rhs(eq)*x/y,
 158    n: ratsimp(x*diff(%g%,x)/(y*diff(%g%,y))),
 159    if not(freeof(x,y,n)) then return(false),
 160    a1: ratsimp(subst(u/x^n,y,%g%)),
 161    a2: integrate(1/(u*(n+a1)),u),
 162    if not(freeof(nounify('integrate),a2)) then return(false),
 163    a3: ratsimp(subst(y*x^n,u,a2)),
 164    method: 'genhom, odeindex: n,
 165    return(x = %c*%e^a3)))$
 166
 167 /*  Chain of solution methods for second order linear homogeneous equations  */
 168
 169 block(local(a1,a2,a3),
 170 hom2(a1,a2,a3):=block([ap,aq,pt],
 171    ap: a2/a1, aq: a3/a1,
 172    if ftest(cc2(ap,aq,y,x)) then return(%q%),
 173    if ftest(exact2(a1,a2,a3)) then return(%q%),
 174    if (pt: pttest(ap)) = false then go(end),
 175    if ftest(euler2(ap,aq)) then return(%q%),
 176    if ftest(bessel2(ap,aq)) then return(%q%),
 177  end,
 178    if ftest(xcc2(ap,aq)) then return(%q%) else return(false)))$
 179
 180 /*  B&DiP, pp. 106-112  */
 181
 182 /* 2nd order linear homogeneous ode with constant coefficients
 183    'diff(y,x,2) + %f%*'diff(y,x) + %g%*y = 0
 184    %f% and %g% independent of x and y */
 185
 186 block(local(%f%,%g%,y,x),
 187 cc2(%f%,%g%,y,x):=block([a,sign,radexpand:'all,alpha,zero,pos,ynew,%k1,%k2],
 188    if not(freeof(x,y,%f%) and freeof(x,y,%g%)) then return(false),
 189    method: 'constcoeff,
 190    a: %f%^2-4*%g%,
 191    if freeof(%i,a) then sign: asksign(a)
 192                    else (radexpand: true, sign: 'pnz),
 193    if sign = zero then
 194      if asksign(%f%^2) = zero then
 195        return(y = %k1 + %k2*x)
 196      else
 197        return(y = %e^(-%f%*x/2) * (%k1 + %k2*x)),
 198    if sign = pos then
 199      return(y = %k1*%e^((-%f%+sqrt(a))*x/2) + %k2*%e^((-%f%-sqrt(a))*x/2)),
 200    a: -a, alpha: x*sqrt(a)/2,
 201    if exponentialize = false then
 202      return(y = %e^(-%f%*x/2) * (%k1*sin(alpha) + %k2*cos(alpha))),
 203    return(y = %e^(-%f%*x/2) * (%k1*exp(%i*alpha) + %k2*exp(-%i*alpha)))))$
 204
 205 /*  B&DiP, pp. 98-99, problem 17  */
 206
 207 block(local(a1,a2,a3),
 208 exact2(a1,a2,a3):=block([b1,%k1,%k2],
 209    if ratsimp(diff(a1,x,2) - diff(a2,x) + a3) = 0
 210      then b1: %e^(-integrate(ratsimp((a2 - diff(a1,x))/a1), x))
 211      else return(false),
 212    method: 'exact,
 213    return(y = %k1*b1*integrate(1/(a1*b1),x) + %k2*b1)))$
 214
 215 /*  B&DiP, pp. 113-114, problem 16  */
 216
 217 block(local(ap,aq),
 218 xcc2(ap,aq):=block([d,b1,z,radexpand:'all],
 219    if aq=0 then return(false),
 220    d: ratsimp((diff(aq,x) + 2*ap*aq)/(2*aq^(3/2))),
 221    if freeof(x,y,d) then b1: cc2(d,1,y,z) else return(false),
 222    method: 'xformtoconstcoeff,
 223    return(subst(integrate(sqrt(aq),x),z,b1))))$
 224 \f
 225 /*  B&DiP, pp. 124-127  */
 226
 227 block(local(soln,%g%),
 228 varp(soln,%g%):=block([y1,y2,y3,y4,wr,heuristic:false,%k1,%k2],
 229    y1: ratsimp(subst([%k1=1,%k2=0],rhs(soln))),
 230    y2: ratsimp(subst([%k1=0,%k2=1],rhs(soln))),
 231    wr: y1*diff(y2,x) - y2*diff(y1,x),
 232    if wr=0 then return(false),
 233    if method='constcoeff and not(freeof('sin,wr)) and not(freeof('cos,wr))
 234      then (heuristic: true, wr: ratsimp(trigreduce(wr))),
 235    y3: ratsimp(y1*%g%/wr),
 236    y4: ratsimp(y2*%g%/wr),
 237    yp: ratsimp(y2*integrate(y3,x) - y1*integrate(y4,x)),
 238    if heuristic=true then yp: ratsimp(trigreduce(yp)),
 239    method: 'variationofparameters,
 240    return(y = rhs(soln) + yp)))$
 241
 242 /*  Methods to reduce second-order equations free of x or y  */
 243
 244 block(local(eq),
 245 reduce(eq):=block([b1,qq],
 246    b1: subst(['diff(y,x,2)=qq, 'diff(y,x)=qq], eq),
 247    if freeof(y,b1) then return(nlx(eq)),
 248    if freeof(x,b1) then return(nly(eq)) else return(false)))$
 249
 250 /*  B&DiP, p. 89, problem 1  */
 251
 252 block(local(eq),
 253 nlx(eq):=block([de,b,a1,v,%k1,%c],
 254    de: subst(['diff(y,x,2)='diff(v,x), 'diff(y,x)=v], eq),
 255    if (b: ode1a(de,v,x)) = false then return(false),
 256    a1: subst([v='diff(y,x),%c=%k1], b),
 257    if ftest(nlxy(a1,'diff(y,x)))
 258      then (method: 'freeofy, return(%q%)) else return(false)))$
 259
 260 /*  B&DiP, p. 89, problem 2  */
 261
 262 block(local(eq),
 263 nly(eq):=block([de,b,a1,yz,v,%c,%k1],
 264    de: subst(['diff(y,x,2)=v*'diff(v,yz), 'diff(y,x)=v, y=yz], eq),
 265    if (b: ode1a(de,v,yz)) = false then return(false),
 266    a1: subst([v='diff(y,x),yz=y,%c=%k1], b),
 267    if ftest(nlxy(a1,'diff(y,x)))
 268      then (method: 'freeofx, return(%q%)) else return(false)))$
 269
 270 block(local(eq,de),
 271 nlxy(eq,de):=block([programmode:true,eq1,%k2,%c],
 272    eq1: solve(eq,de),
 273    eq1: maplist(lambda([zz], if ftest(ode1a(zz,y,x))
 274                                then subst(%k2,%c,%q%) else false),
 275                 eq1),
 276    if length(eq1)=1 then return(first(eq1)) else return(eq1)))$
 277 \f
 278 /*    This is a start on a series of programs to recognize and
 279     solve certain special classes of differential equations.
 280     In particular, to start with, is the Euler, or equidimensional,
 281     equation:  x^2*y'' + axy' + by = 0.  Actually, the form we
 282     will investigate is:  y'' + ay'/x + by/x^2 = 0.
 283       PTTEST analyzes the y' term for a coefficient of the form
 284     a/(x-pt), since we must assume that the equation may be
 285     translated.  */
 286
 287 block(local(a),
 288 pttest(a):=block([a1,a2,a3],
 289    if (a1: ratsimp(a)) = 0 then return(false),
 290    a1: expand(1/a1),
 291    if (a2: coeff(a1,x,1)) = 0 then return(false),
 292    if not(freeof(x,a2)) then return(false),
 293    a3: coeff(a1,x,0),
 294    if not(a1 = a2*x + a3) then return(false) else return(-a3/a2)))$
 295
 296 block(local(a,b),
 297 euler2(a,b):=block([dc,rp,ip,alpha,beta,sign,radexpand:false,%k1,%k2,pos,zero],
 298    if not(freeof(x,y,beta: ratsimp(b*(x-pt)^2))) then return(false),
 299    method: 'euler, alpha: a*(x-pt),
 300    dc: ratsimp((alpha-1)^2 - 4*beta),
 301    rp: ratsimp(-(alpha-1)/2),
 302    sign: asksign(dc),
 303    if sign = zero then return(y = (x-pt)^rp * (%k1 + %k2*log(x-pt))),
 304    if sign = pos then
 305      (ip: sqrt(dc)/2, return(y = %k1*(x-pt)^(rp+ip) + %k2*(x-pt)^(rp-ip))),
 306    dc: -dc, ip: sqrt(dc)/2,
 307    return(y = (x-pt)^rp * (%k1*sin(ip*log(x-pt)) + %k2*cos(ip*log(x-pt))))))$
 308
 309 block(local(a,b),
 310 bessel2(a,b):=block([nu,b1,intp,radexpand:'all,%k1,%k2],
 311    if not(freeof(x,y,b1: ratsimp((1-b)*(x-pt)^2))) then return(false),
 312    if ratsimp(a*(x-pt)) # 1 then return(false),
 313    nu: sqrt(b1), method: 'bessel,
 314    if nu = 1/2 then return(y = (%k1*sin(x-pt) + %k2*cos(x-pt))/sqrt(x-pt)),
 315    intp:askinteger(nu),
 316    if intp = 'yes then return(y = %k1*bessel_j(nu,x-pt) + %k2*bessel_y(nu,x-pt)),
 317    return(y = %k1*bessel_j(nu,x-pt) + %k2*bessel_j(-nu,x-pt))))$
 318
 319 block(local(soln,xc,yc),
 320 ic1(soln,xc,yc):=
 321    block([%c],
 322    (noteqn(xc), noteqn(yc), boundtest('%c,%c),
 323    ratsimp(subst(['%c=rhs(solve1(at(soln,[xc,yc]),%c))],soln)))))$
 324
 325 block(local(soln,xa,ya,xb,yb),
 326 bc2(soln,xa,ya,xb,yb):=
 327    block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2],
 328       noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb),
 329       boundtest('%k1,%k1), boundtest('%k2,%k2),
 330       temp: maplist(lambda([zz], subst(zz,soln)),
 331                     solve([subst([xa,ya],soln),
 332                            subst([xb,yb],soln)],
 333                           [%k1,%k2])),
 334       if length(temp)=1 then return(first(temp)) else return(temp)))$
 335
 336 block(local(soln,xa,ya,dya),
 337 ic2(soln,xa,ya,dya):=
 338    block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1,yas],
 339       noteqn(xa), noteqn(ya), noteqn(dya),
 340       boundtest('%k1,%k1), boundtest('%k2,%k2),
 341       yas:lhs(ya),
 342       /* We need a dependency for the algorithm. Add it, if it is not present */
 343       if diff(yas,lhs(xa)) = 0 then depends(yas,lhs(xa)),
 344       if listp(soln) then
 345          /* A list of solutions. Map over the solutions. */
 346          temp: map('lhs,soln)-map('rhs, soln)
 347       else
 348          temp: lhs(soln) - rhs(soln),
 349       temp: maplist(lambda([zz], subst(zz,soln)),
 350                     solve(flatten([subst([xa,ya], soln),
 351                                    subst([dya,xa,ya], diff(soln,lhs(xa)))]),
 352                           [%k1,%k2])),
 353       if length(temp)=1 then return(first(temp)) else return(temp)))$
 354
 355 block(local(x),
 356 noteqn(x):=if atom(x) or not inpart(x,0)="="
 357              then (disp(x), disp("not an equation"), error()))$
 358
 359 block(local(x,y),
 360 boundtest(x,y):=
 361    if x#y then (disp(x), disp("must not be bound"), error()))$
 362
 363 block(local(msg,eq),
 364 failure(msg,eq):=
 365    block([ynew],   (if not status(feature,"ode")
 366       then (ldisp(subst(yold,ynew,eq)), disp(msg)),
 367     false)))$
 368
 369 msg1:   "not a proper differential equation"$
 370 msg2:   "first order equation not linear in y'"$
 371