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