share/algebra/recur.mac

   1 define_variable(g,'g,any)$
   2
   3 /*this block checks for a polynomial in n*/
   4 polyp(g,n):=block([d,f,c],
   5      g:ratexpand(g), if freeof(n,g) then return(true),
   6      d:hipow(g,n), f:true,
   7      for i:d step -1 thru 0 do
   8      (c:coeff(g,n,i), if not(freeof(n,c)) then f:false,
   9       g:ratexpand(g-c*n^i)),
  10      return(is(g=0 and f)))$
  11
  12
  13
  14 /*this block checks for a constant to a polynomial power*/
  15 polyinn(x,n):=block([b,e],
  16      if inpart(x,0)="*" then
  17      return(polyinn(inpart(g,1),n) and polyinn(inpart(g,2),n)),
  18      if inpart(x,0)#"^" then return(false),
  19      b:inpart(x,1),
  20      e:inpart(x,2),
  21      if not freeof(n,b) then return(false),
  22      return(polyp(e,n)))$
  23
  24
  25
  26 /*this block implements the characteristic equation method*/
  27 char(e,g,u,n,k,iv):=block([gensol,homsol,parsol,los,multiplicities,
  28 h,v,l,ss,dispflag],
  29 local(a,aa,b,r,m),
  30 dispflag:false,
  31   for i:0 thru k do
  32     aa[i]:coeff(e,u(n+k-i)),
  33     h:0,
  34   for i:0 thru k do
  35     h:h+aa[i]*u(n+k-i),
  36   if h#e then return("erroneous input"),
  37
  38   for i:0 thru k do
  39     h:subst(u^(k-i),u(n+k-i),h),
  40
  41   multiplicities:true,
  42   los:solve(h,u),
  43   for i:1 thru length(los) do
  44     (r[i]:los[i], r[i]:rhs(ev(r[i])),
  45      m[i]:multiplicities[i]),
  46   homsol:
  47     sum(sum(a[i,j]*n^(m[i]-j),j,1,m[i])*r[i]^n,i,1,length(los)),
  48
  49
  50   if g=0 then
  51      (v:[ ],
  52       for i:1 thru length(los) do
  53       for j:1 thru m[i] do v:cons(a[i,j],v),
  54          l:[ ],
  55       for q:0 thru k-1 do l:cons(subst(q,n,homsol)=u(q),l),
  56       ss:ev(solve(l,v),iv),
  57       return(u(n)=(ev(homsol,ss))))
  58
  59   else if polyp(g,n) = true then
  60      (g:ratexpand(g), parsol:sum(b[j]*n^j,j,0,hipow(g,n)),
  61       for j:0 thru k do
  62        (l:0, v:e,
  63          for i:0 thru k do
  64          (l:ratexpand(subst(n+k-i,n,b[j]*n^j)),
  65           v:ratexpand(subst(l,u(n+k-i),v))),
  66           v:ratsimp(v),
  67            if v#0 then return(v) else parsol:n*parsol),
  68          v:e,
  69       for i:0 thru k do (l:ratexpand(subst(n+k-i,n,parsol)),
  70          v:ratexpand(subst(l,u(n+k-i),v))),
  71          l:[ ],
  72       for i:0 thru hipow(parsol,n) do
  73          l:cons(coeff(v=g,n,i),l),
  74          v:[ ],
  75       for j:0 thru hipow(parsol,n) do
  76          v:cons(b[j],v),
  77       ss:solve(l,v),
  78       parsol:ev(parsol,ss))
  79
  80   else if polyinn(g,n) = true then
  81      (parsol:b1*g,
  82       for j:0 thru k do
  83        (l:0, v:e,
  84          for i:0 thru k do
  85          (l:subst(n+k-i,n,parsol), v:subst(l,u(n+k-i),v)),
  86           v:ratsimp(v),
  87           if v#0 then return(v) else parsol:n*parsol),
  88       ss:solve(v=g,b1),
  89       parsol:ev(parsol,ss))
  90
  91   else if inpart(g,0)=sin or inpart(g,0) = cos then
  92      (parsol:b[1]*sin(inpart(g,1)) + b[2]*cos(inpart(g,1)),
  93       for j:0 thru k do
  94        (l:0, v:e,
  95          for i:0 thru k do
  96          (l:expand(subst(n+k-i,n,parsol)),
  97           v:expand(subst(l,u(n+k-i),v))),
  98       v:trigexpand(v),
  99       if v#0 then return(v) else parsol:n*parsol),
 100          v:e,
 101       for i:0 thru k do(l:expand(subst(n+k-i,n,parsol)),
 102          v:expand(subst(l,u(n+k-i),v))),
 103       v:trigexpand(v),
 104          l:[ ],
 105       lt:[sin(inpart(g,1)),cos(inpart(g,1))],
 106       for jj:1 thru 2 do
 107          l:cons(coeff(v=g,lt[jj]),l),
 108          v:[ ],
 109       for j:1 thru 2 do
 110          v:cons(b[j],v),
 111       ss:solve(l,v),
 112       parsol:ev(parsol,ss))
 113
 114
 115   else return("can't be solved in closed form by program"),
 116
 117   gensol:homsol + parsol,
 118      v:[ ],
 119   for i:1 thru length(los) do
 120   for j:1 thru m[i] do v:cons(a[i,j],v),
 121      l:[ ],
 122   for q:0 thru k-1 do
 123   l:cons(subst(q,n,gensol)=u(q),l),
 124   ss:ev(solve(l,v),iv),
 125   return(u(n)=(ev(gensol,ss))))$
 126
 127
 128
 129
 130 /*this block implements the generating function method*/
 131 genf(e,g,u,n,k,iv):=block([multiplicities,l,v,ss,vv,los,
 132 nr,f,sol,p,dispflag],
 133 local(a,aa,b),
 134 dispflag:false,
 135   for i:0 thru k do
 136     aa[i]:coeff(e,u(n+k-i)),
 137     h:0,
 138   for i:0 thru k do
 139     h:h+aa[i]*u(n+k-i),
 140   if h#e then return("erroneous input"),
 141
 142      l:e,
 143   for i:0 thru k do
 144     l:subst((f-sum(u(j)*x^j,j,0,k-i-1))*x^i,u(n+k-i),l),
 145
 146   if g=0 then
 147      (s:solve(l,f),
 148       f:ev(f,s))
 149
 150   else if polyp(g,n) = true then
 151      (g:ratexpand(g),
 152       v:subst(x^k/(1-x)*coeff(g,n,0),coeff(g,n,0),g),
 153       vv:ratsimp(diff(1/(1-x),x)),
 154       for i:1 thru hipow(g,n) do
 155         (v:subst(x^k*x*vv*coeff(g,n,i),coeff(g,n,i)*n^i,v),
 156          vv:ratsimp(diff(x*vv,x))),
 157       v:ratsimp(v),
 158       ss:solve(l=v,f),
 159       f:ev(f,ss))
 160
 161   else if polyinn(g,n) = true and hipow(inpart(g,2),n) < 2 then
 162      (g1:(x^k)*(inpart(g,1)^coeff(inpart(g,2),n,0)),
 163       g2:1 - x*(inpart(g,1)^coeff(inpart(g,2),n,1)),
 164       v:ratsimp(g1/g2),
 165       ss:solve(l=v,f),
 166       f:ev(f,ss))
 167
 168   else return("can't be solved in closed form by program"),
 169
 170
 171   multiplicities:true,
 172   los:solve(newrat(f),x),
 173   for i:1 thru length(los) do
 174      (r[i]:los[i], r[i]:rhs(ev(r[i])),
 175       m[i]:multiplicities[i]),
 176
 177      v:[ ],
 178   b:product((1-r[i]*x)^m[i],i,1,length(los)),
 179   for i:1 thru length(los) do
 180   for j:1 thru m[i] do
 181      (p[i,j]:b*a[i,j]/((1-r[i]*x)^j),  v:cons(a[i,j],v)),
 182   p:sum(sum(p[i,j],j,1,m[i]),i,1,length(los)),
 183
 184      l:[ ],
 185   nf:ratexpand(num(f)/abs(coeff(denom(f),x,0))), p:ratexpand(p),
 186   for i:0 thru hipow(ratexpand(b),x)-1  do
 187      l:cons(coeff(nf=p,x,i),l),
 188   sss:ev(solve(l,v),iv),
 189
 190   sol:sum(sum(a[i,j]*coeff(denom(f),x,0)/abs(coeff(denom(f),x,0))*
 191      binomial(j+n-1,n)*r[i]^n,j,1,m[i]),i,1,length(los)),
 192
 193   return(u(n)=(ev(sol,sss))))$
 194
 195
 196 /*this block finds the new polynomial associated to f*/
 197 newrat(f):=block([hd,cp,dp],
 198      hd:hipow(denom(f),x),
 199      cp:coeff(denom(f),x,hd),
 200      dp:sum((coeff(denom(f),x,i))/cp*x^i,i,0,hd),
 201      return(sum(coeff(dp,x,hd-i)*x^i,i,0,hd)))$
 202
 203
 204
 205
 206
 207 /*this block implements the variable coefficient method*/
 208 varc1(e,g,u,n,k,iv):=block([v,vv,eq,y,cauchysum,finsol,sersol,dispflag],
 209 local(a,b),dispflag:false,
 210   for i:0 thru k do
 211     (a[i]:coeff(e,u(n+i)),
 212      a[i]:ratexpand(a[i]),
 213      if polyp(a[i],n)=false then return("can't do it")),
 214   if k=2 and (b:besselcheck(e,k) # false) then return(b),
 215   v:ratexpand(e),
 216   for i:k step -1 thru 0 do
 217      for j:hipow(a[i],n) step -1 thru 0 do
 218        (v:ratsubst(x^j*'diff(y,x,i+j),n^j*u(n+i),v),
 219         v:ratexpand(v)),
 220   v:ratsubst(y,'diff(y,x,0),v),
 221   v:ratexpand(v),
 222   if polyp(g,n) = true then
 223     (g:ratexpand(g), vv:g,
 224      for i:0 thru hipow(g,n) do
 225         vv:subst(x^i,n^i,vv),
 226         vv:%e^x*vv)
 227   else return("can't do it"),
 228
 229      eq:v-vv,
 230   dependencies(y(x)),
 231   sol:ode2(eq=0,y,x),
 232         if k=1 then finsol:ic1(sol,x=0,y=ev(u(0),iv))
 233         else if k=2 then finsol:ic2(sol,x=0,y=ev(u(0),iv),'diff(y,x)=ev(u(1),iv))
 234   else return("o.d.e. can't be solved at present by macsyma"),
 235      cauchysum:true,
 236   sersol:powerseries(rhs(finsol),x,0), sersol:expand(sersol),
 237   b:inpart(sersol,2),
 238   b:ev(b,x=1),
 239   if atom(b)=false then b:substpart(n,b,4),
 240   return(u(n)=((n!)*b)))$
 241
 242
 243 /*this block checks for a bessel recurrence relation*/
 244
 245 besselcheck(e,k):=block([a,ans],
 246 local(a),
 247      for i:0 thru k do
 248      (a[i]:coeff(e,u(n+i)),
 249       a[i]:ratexpand(a[i])),
 250      if not(integerp(a[0])) then return(false),
 251      if not(integerp(ev(a[1],n=0))) then return(false),
 252      if not(hipow(a[1],n)=1) then return(false),
 253      if not(integerp(coeff(a[1],n,1))) then return(false),
 254      if not(a[2]=1) then return(false),
 255      ans:"a linear combination of bessel functions",
 256      /*exact details are of no significance,since we are merely
 257        demonstrating the feasibility of this approach*/
 258      return(ans))$
 259
 260
 261
 262
 263
 264 /*this block implements the first order method*/
 265 varc2(e,g,u,n,k,iv):=block([h,p,v,c,sol,sj,si,pn,hi,ap],
 266      p:(-1)*coeff(e,u(n))/coeff(e,u(n+1)),
 267      v:g/coeff(e,u(n+1)),
 268      sj:subst(j,n,p),
 269      si:subst(i,n,p),
 270      pn:product(si,i,1,n-1),
 271      hi:subst(i,n,v)/product(sj,j,1,i-1),
 272      v1:sum(hi,i,0,n),
 273      ap:ev(u(0)-subst(0,n,v),iv),
 274      return(u(n)=ap*pn+pn*v1))$