/* -*- Mode: MACSYMA; Package: MAXIMA -*- */

/*  (c) Copyright 1981 Massachusetts Institute of Technology  */

/*  The Ordinary Differential Equation Solver.
    This package consists primarily of a set of routines taken from Moses'
    thesis and Boyce & DiPrima for solving O.D.E.s of 1st and 2nd order. 
    The top-level routines are ODE2, IC1, IC2, and BC2.  */

/* The following declares help the translated code.  We also used the
   convention  of %f% being the special rather than declaring f special
   since other programs might use f */

eval_when(translate,
declare_translated(boundtest,noteqn,nlxy,nly,nlx,xcc2,bessel2,euler2,
   pttest,exact2,cc2,genhom,solvebernoulli,solvehom,integfactor,exact,
   separable,solvelnr,solve1,linear2,reduce,hom2,pr2,varp,desimp,failure,
   ode1a,ftest,ode2a));

/* Remove some of the specials like x and y at the end */
declare(pt,special,yp,special,yold,special,%q%,special,x,special,y,special,
   method,special,%f%,special,
    %g%,special,msg1,special,
    msg2,special,intfactor,special,odeindex,special,singsolve,special);

ode2(eq,yold,x):=block([derivsubst:false,ynew],
   subst(yold,ynew,ode2a(subst(ynew,yold,eq),ynew,x)))$

ode2a(eq,y,x):=block([de,a1,a2,a3,a4,%q%,msg1],
   intfactor: false, method: 'none,
   if freeof('diff(y,x,2),eq)
     then if ftest(ode1a(eq,y,x)) then return(%q%) else return(false),
   if derivdegree(de: desimp(lhs(eq)-rhs(eq)),y,x) # 2
     then return(failure(msg1,eq)),
   a1: coeff(de,'diff(y,x,2)),
   a2: coeff(de,'diff(y,x)),
   a3: coeff(de,y),
   a4: expand(de - a1*'diff(y,x,2) - a2*'diff(y,x) - a3*y),
   if pr2(a1) and pr2(a2) and pr2(a3) and pr2(a4) and
      ftest(hom2(a1,a2,a3))
     then if a4=0 then return(%q%) else return(varp(%q%,-a4/a1)),
   if ftest(reduce(de)) then return(%q%) else return(false))$

ode1a(eq,y,x):=block([de,des], /* f, g, %q% */
   if derivdegree(de: expand(lhs(eq)-rhs(eq)),y,x) # 1
     then return(failure(msg1,eq)),
   if linear2(de,'diff(y,x)) = false then return(failure(msg2,eq)),
   des: desimp(de),
   de: solve1(des,'diff(y,x)),
   if ftest(solvelnr(de)) then return(%q%),
   if ftest(separable(de)) then return(%q%),
   if ftest(integfactor(%g%,%f%)) then return(exact(%q%*%g%,%q%*%f%)),
			       /* linear2 binds %f% and %g% */
   if linear2(des,'diff(y,x)) = false then return(failure(msg2,eq)),
   if ftest(integfactor(%g%,%f%)) then return(exact(%q%*%g%,%q%*%f%)),
   if ftest(solvehom(de)) then return(%q%),
   if ftest(solvebernoulli(de)) then return(%q%),
   if ftest(genhom(de)) then return(%q%) else return(false))$

desimp(eq):=block([inflag:true],
   eq: factor(eq),
   if atom(eq) or not(inpart(eq,0)="*") then return(expand(eq)),
   eq: map(lambda([u], if freeof(nounify('diff),u) then 1 else u), eq),
   return(expand(eq)))$

pr2(%f%):=freeof(y,'diff(y,x),'diff(y,x,2),%f%)$

ftest(call):=is((%q%: call) # false)$

solve1(eq,y):=block([programmode:true],first(solve(eq,y)))$

/* linear2 tests for the form fx+%g% */

linear2(expr,x):=block([],
   %f%: ratcoef(expr,x),
   if not(freeof(x,%f%)) then return(false),
   %g%: ratsimp(expr - %f%*x),
   return(freeof(x,%g%)))$

/*  variables used to denote constants: %C, %K1, %K2.
    METHOD denotes the method of solution.
    INTFACTOR denotes the integrating factor.
    ODEINDEX denotes the index for Bernoulli's method or for the genhom method.
    YP denotes the particular solution for the variation of parameters technique.  */

/*  B&DiP, pp. 13-14  */

solvelnr(eq):=block([%f%,%g%,w,%c],
   if linear2(rhs(eq),y) = false then return(false),
   w: %e^(integrate(%f%,x)),
   method: 'linear,
   return(y=w*(integrate(%g%/w,x)+%c)))$

/*  B&DiP, pp. 29-34  */

separable(eq):=block([xpart:[],ypart:[],flag:false,inflag:true,%c],
   eq: factor(rhs(eq)),
   if atom(eq) or not(inpart(eq,0)="*") then eq: [eq],
   for u in eq do
      if freeof(x,u) then ypart: cons(u,ypart) else
      if freeof(y,u) then xpart: cons(u,xpart) else return(flag: true),
   if flag = true then return(false),
   if xpart = [] then xpart: 1 else xpart: apply("*",xpart),
   if ypart = [] then ypart: 1 else ypart: apply("*",ypart),
   method: 'separable,
   return(ratsimp(integrate(1/ypart,y)) = ratsimp(integrate(xpart,x)) + %c))$

/*  B&DiP, pp. 34-41  */

integfactor(m,n):=block([b1,b2,dmdx,dmdy,dndx,dndy,dd,%e_to_numlog:true],
   dmdy: ratsimp(diff(m,y)),  dndx: ratsimp(diff(n,x)),
   if (dd: dmdy-dndx) = 0 then return(1),
   dmdx: ratsimp(diff(m,x)),  dndy: ratsimp(diff(n,y)),
   if dmdx-dndy=0 and dmdy+dndx=0 then return(1/(m^2 + n^2)),
   if freeof(y, (b1: ratsimp(dd/n))) then return(%e^(integrate(b1,x))),
   if freeof(x, (b2: ratsimp(dd/m)))
     then return(%e^(integrate(-b2,y))) else return(false))$

exact(m,n):=block([a,ynew,%c],
   intfactor: subst(yold,ynew,%q%),
   a: integrate(ratsimp(m),x),
   method: 'exact,
   return(ratsimp(a + integrate(ratsimp(n-diff(a,y)),y)) = %c))$

/*  B&DiP, pp. 43-44  */

solvehom(eq):=block([qq,a1,a2,%c],
   a1: ratsimp(subst(x*qq,y,rhs(eq))),
   if not(freeof(x,a1)) then return(false),
   a2: ratsimp(subst(y/x,qq,integrate(1/(a1-qq),qq))),
   method: 'homogeneous,
   return(%c*x = %e^a2))$

/*  B&DiP, p. 21, problem 15  */

solvebernoulli(eq):=block([a1,a2,n,%c],
   a1: coeff(eq: expand(rhs(eq)),y,1),
   if not(freeof(y,a1)) then return(false),
   n: hipow(ratsimp(eq-a1*y),y),
   a2: coeff(eq,y,n),
   if not(freeof(y,a2)) or not(freeof(x,y,n)) or n=0
		        or not(eq = expand(a1*y + a2*y^n))
     then return(false),
   a1: integrate(a1,x),
   method: 'bernoulli, odeindex: n,
   return(y = %e^a1 * ((1-n)*integrate(a2*%e^((n-1)*a1),x) + %c) ^ (1/(1-n))))$

/*  Generalized homogeneous equation:  y' = y/x * H(yx^n)
	Reference:  Moses' thesis.  */

genhom(eq):=block([%g%,u,n,a1,a2,a3,%c],
   %g%: rhs(eq)*x/y,
   n: ratsimp(x*diff(%g%,x)/(y*diff(%g%,y))),
   if not(freeof(x,y,n)) then return(false),
   a1: ratsimp(subst(u/x^n,y,%g%)),
   a2: integrate(1/(u*(n+a1)),u),
   if not(freeof(nounify('integrate),a2)) then return(false),
   a3: ratsimp(subst(y*x^n,u,a2)),
   method: 'genhom, odeindex: n,
   return(x = %c*%e^a3))$

/*  Chain of solution methods for second order linear homogeneous equations  */

hom2(a1,a2,a3):=block([ap,aq,pt],
   ap: a2/a1, aq: a3/a1, 
   if ftest(cc2(ap,aq,y,x)) then return(%q%),
   if ftest(exact2(a1,a2,a3)) then return(%q%),
   if (pt: pttest(ap)) = false then go(end),
   if ftest(euler2(ap,aq)) then return(%q%),
   if ftest(bessel2(ap,aq)) then return(%q%),
 end, 
   if ftest(xcc2(ap,aq)) then return(%q%) else return(false))$

/*  B&DiP, pp. 106-112  */

cc2(%f%,%g%,y,x):=block([a,sign,radexpand:'all,alpha,zero,pos,ynew,%k1,%k2],
   if not(freeof(x,y,%f%) and freeof(x,y,%g%)) then return(false),
   method: 'constcoeff,
   a: %f%^2-4*%g%,
   if freeof(%i,a) then sign: asksign(a)
		   else (radexpand: true, sign: 'pnz),
   if sign = zero then return(y = %e^(-%f%*x/2) * (%k1 + %k2*x)),
   if sign = pos then
     return(y = %k1*%e^((-%f%+sqrt(a))*x/2) + %k2*%e^((-%f%-sqrt(a))*x/2)),
   a: -a, alpha: x*sqrt(a)/2,
   if exponentialize = false then
     return(y = %e^(-%f%*x/2) * (%k1*sin(alpha) + %k2*cos(alpha))),
   return(y = %e^(-%f%*x/2) * (%k1*exp(%i*alpha) + %k2*exp(-%i*alpha))))$ 

/*  B&DiP, pp. 98-99, problem 17  */

exact2(a1,a2,a3):=block([b1,%k1,%k2],
   if ratsimp(diff(a1,x,2) - diff(a2,x) + a3) = 0
     then b1: %e^(-integrate(ratsimp((a2 - diff(a1,x))/a1), x))
     else return(false),
   method: 'exact,
   return(y = %k1*b1*integrate(1/(a1*b1),x) + %k2*b1))$

/*  B&DiP, pp. 113-114, problem 16  */

xcc2(ap,aq):=block([d,b1,z,radexpand:'all],
   if aq=0 then return(false),
   d: ratsimp((diff(aq,x) + 2*ap*aq)/(2*aq^(3/2))),
   if freeof(x,y,d) then b1: cc2(d,1,y,z) else return(false),
   method: 'xformtoconstcoeff,
   return(subst(integrate(sqrt(aq),x),z,b1)))$

/*  B&DiP, pp. 124-127  */

varp(soln,%g%):=block([y1,y2,y3,y4,wr,heuristic:false,%k1,%k2],
   y1: ratsimp(subst([%k1=1,%k2=0],rhs(soln))),
   y2: ratsimp(subst([%k1=0,%k2=1],rhs(soln))),
   wr: y1*diff(y2,x) - y2*diff(y1,x),
   if wr=0 then return(false),
   if method='constcoeff and not(freeof('sin,wr)) and not(freeof('cos,wr))
     then (heuristic: true, wr: ratsimp(trigreduce(wr))),
   y3: ratsimp(y1*%g%/wr),
   y4: ratsimp(y2*%g%/wr),
   yp: ratsimp(y2*integrate(y3,x) - y1*integrate(y4,x)),
   if heuristic=true then yp: ratsimp(trigreduce(yp)),
   method: 'variationofparameters,
   return(y = rhs(soln) + yp))$

/*  Methods to reduce second-order equations free of x or y  */

reduce(eq):=block([b1,qq],
   b1: subst(['diff(y,x,2)=qq, 'diff(y,x)=qq], eq),
   if freeof(y,b1) then return(nlx(eq)),
   if freeof(x,b1) then return(nly(eq)) else return(false))$

/*  B&DiP, p. 89, problem 1  */

nlx(eq):=block([de,b,a1,v,%k1,%c],
   de: subst(['diff(y,x,2)='diff(v,x), 'diff(y,x)=v], eq),
   if (b: ode1a(de,v,x)) = false then return(false),
   a1: subst([v='diff(y,x),%c=%k1], b),
   if ftest(nlxy(a1,'diff(y,x)))
     then (method: 'freeofy, return(%q%)) else return(false))$

/*  B&DiP, p. 89, problem 2  */

nly(eq):=block([de,b,a1,yz,v,%c,%k1],
   de: subst(['diff(y,x,2)=v*'diff(v,yz), 'diff(y,x)=v, y=yz], eq),
   if (b: ode1a(de,v,yz)) = false then return(false),
   a1: subst([v='diff(y,x),yz=y,%c=%k1], b),
   if ftest(nlxy(a1,'diff(y,x)))
     then (method: 'freeofx, return(%q%)) else return(false))$

nlxy(eq,de):=block([programmode:true,eq1,%k2,%c],
   eq1: solve(eq,de),
   eq1: maplist(lambda([zz], if ftest(ode1a(zz,y,x))
			       then subst(%k2,%c,%q%) else false),
		eq1),
   if length(eq1)=1 then return(first(eq1)) else return(eq1))$

/*    This is a start on a series of programs to recognize and 
    solve certain special classes of differential equations. 
    In particular, to start with, is the Euler, or equidimensional, 
    equation:  x^2*y'' + axy' + by = 0.  Actually, the form we 
    will investigate is:  y'' + ay'/x + by/x^2 = 0.
      PTTEST analyzes the y' term for a coefficient of the form 
    a/(x-pt), since we must assume that the equation may be 
    translated.  */

pttest(a):=block([a1,a2,a3],
   if (a1: ratsimp(a)) = 0 then return(false),
   a1: expand(1/a1),
   if (a2: coeff(a1,x,1)) = 0 then return(false),
   if not(freeof(x,a2)) then return(false),
   a3: coeff(a1,x,0),
   if not(a1 = a2*x + a3) then return(false) else return(-a3/a2))$

euler2(a,b):=block([dc,rp,ip,alpha,beta,sign,radexpand:false,%k1,%k2,pos,zero],
   if not(freeof(x,y,beta: ratsimp(b*(x-pt)^2))) then return(false),
   method: 'euler, alpha: a*(x-pt),
   dc: ratsimp((alpha-1)^2 - 4*beta),
   rp: ratsimp(-(alpha-1)/2),
   sign: asksign(dc),
   if sign = zero then return(y = (x-pt)^rp * (%k1 + %k2*log(x-pt))),
   if sign = pos then 
     (ip: sqrt(dc)/2, return(y = %k1*(x-pt)^(rp+ip) + %k2*(x-pt)^(rp-ip))),
   dc: -dc, ip: sqrt(dc)/2,
   return(y = (x-pt)^rp * (%k1*sin(ip*log(x-pt)) + %k2*cos(ip*log(x-pt)))))$

bessel2(a,b):=block([nu,b1,intp,radexpand:'all,%k1,%k2],
   if not(freeof(x,y,b1: ratsimp((1-b)*(x-pt)^2))) then return(false),
   if ratsimp(a*(x-pt)) # 1 then return(false),
   nu: sqrt(b1), method: 'bessel,
   if nu = 1/2 then return(y = (%k1*sin(x-pt) + %k2*cos(x-pt))/sqrt(x-pt)),
   if featurep(nu,'integer) then intp: 'y else if numberp(nu) then intp: 'n,
 loop,
   if not(intp='y or intp='n)
     then (intp: readonly("is", nu, "an integer?  type y or n."), go(loop)),
   if intp = 'y then return(y = %k1*bessel_j(nu,x-pt) + %k2*bessel_y(nu,x-pt)),
   return(y = %k1*bessel_j(nu,x-pt) + %k2*bessel_j(-nu,x-pt)))$

ic1(soln,xc,yc):=
   block([%c],
   (noteqn(xc), noteqn(yc), boundtest('%c,%c),
    ratsimp(subst(['%c=rhs(solve1(subst([xc,yc],soln),%c))],soln))))$

bc2(soln,xa,ya,xb,yb):=
   block([programmode:true,backsubst:true,singsolve:true,temp,%k1,%k2],
      noteqn(xa), noteqn(ya), noteqn(xb), noteqn(yb),
      boundtest('%k1,%k1), boundtest('%k2,%k2),
      temp: maplist(lambda([zz], subst(zz,soln)),
		    solve([subst([xa,ya],soln),
			   subst([xb,yb],soln)],
			  [%k1,%k2])),
      if length(temp)=1 then return(first(temp)) else return(temp))$

ic2(soln,xa,ya,dya):=
   block([programmode:true,backsubst:true,singsolve:true,temp,%k2,%k1],
      noteqn(xa), noteqn(ya), noteqn(dya),
      boundtest('%k1,%k1), boundtest('%k2,%k2),
      temp: lhs(soln) - rhs(soln),
      temp: maplist(lambda([zz], subst(zz,soln)),
		    solve([subst([xa,ya],soln), subst([dya,xa],
			     lhs(dya)=-subst(0,lhs(dya),diff(temp,lhs(xa)))
				       /diff(temp,lhs(ya)))],
			  [%k1,%k2])),
      if length(temp)=1 then return(first(temp)) else return(temp))$

noteqn(x):=if atom(x) or not inpart(x,0)="="
	     then (disp(x), disp("not an equation"), error())$

boundtest(x,y):=
   if x#y then (disp(x), disp("must not be bound"), error())$

failure(msg,eq):=
   block([ynew],   (if not status(feature,"ode")
      then (ldisp(subst(yold,ynew,eq)), disp(msg)),
    false))$

eval_when([translate,load,compile],
       remove(x,special,y,special))$
 
msg1:	"not a proper differential equation"$
msg2:	"first order equation not linear in y'"$