share/contrib/rand/reduct3.mac

   1 /* Filename reduct3.mac
   2
   3    ***************************************************************
   4    *                                                             *
   5    *                     <package name>                          *
   6    *                <functionality description>                  *
   7    *                                                             *
   8    *          from: Perturbation Methods, Bifurcation            *
   9    *                Theory and Computer Algebra.                 *
  10    *           by Rand & Armbruster (Springer 1987)              *
  11    *                Programmed by Richard Rand                   *
  12    *      These files are released to the public domain          *
  13    *                                                             *
  14    ***************************************************************
  15 */
  16 /* program number 12: a liapunov-schmidt reduction for steady state
  17    bifurcations in systems of partial differential equations
  18    depending on one independent space variable. see page 189 in
  19    "perturbation methods, bifurcation theory and computer algebra".
  20 */
  21
  22
  23
  24
  25 /*this file contains reduction3(), a function to perform a liapunov-
  26 schmidt reduction for steady state bifurcations from n coupled
  27 partial differential equations on one spatial dimension.
  28 the following additional functions are supplied:
  29   setup() allows the problem to be entered.
  30   g_poly(i,j) calculates the coefficient of amp^i lam^j in the bifurcation
  31   equation g(amp,lam).
  32   w_poly(i,j) calculates the coefficient of amp^i lam^j in w(x;amp,lam).
  33   solve_ode(exp) solves certain ordinary differential equations via a fourier
  34   mode ansatz.
  35   feedw(exp) ensures that <afun,w> = 0 .
  36   find_trig(exp) identifies fourier modes.
  37   setify(list) transforms a list into a set.
  38   g_eq() assembles the bifurcation equation.
  39   vfun(list,value) creates the substitution list
  40           [list[1] = value, list[2] = value, ...]
  41   diffnull(i,j) sets the derivative diff(w,amp,i,lam,j) to zero.*/
  42
  43 reduction3():=block(
  44        setup(),
  45        order:read("to which order do you want to calculate"),
  46        for i:2 thru order-1 do  w_poly(i,0),
  47        for i:1 thru order-2 do  w_poly(i,1),
  48        for i:1 thru order do g_poly(i,0),
  49        for i:1 thru order-1 do g_poly(i,1),
  50        g_eq())$
  51
  52
  53 setup():=( /* this function performs the input for the liapunov-schmidt
  54               reduction*/
  55 assume_pos:true,
  56 ls_list:[],
  57 n:read("enter the number of differential equations"),
  58 y:read("enter the dependent variables as a list"),
  59 xvar:read("enter the spatial coordinate"),
  60 alpha:read("enter the bifurcation parameter"),
  61 cal:read("enter the critical bifurcation value"),
  62 print("we define lam = ",alpha - cal),
  63 cfun:read("enter the critical eigenfunction as a list"),
  64 afun:read("enter the adjoint critical eigenfunction as a list"),
  65 kill(w),
  66 w:makelist(concat(w,i),i,1,n),
  67 zwlist:makelist(concat(zw,i),i,1,n),
  68 nullist:makelist(0,i,1,n),
  69 depends(append(zwlist,w,y),cons(xvar,[amp,lam])),
  70 eqs:makelist(read("enter the differential equation number",i),i,1,n),
  71 eqlam:ev(eqs,ev(alpha) = lam + cal),
  72 print(eqlam),
  73 len:read("what is the length of the space interval"),
  74 wnull:vfun(w,0),
  75 sub:maplist("=",y,amp*cfun+w),
  76 database:append(difnull(1,0),difnull(0,1)),
  77 zwnull:vfun(zwlist,0),
  78 norm:integrate(afun.cfun,xvar,0,len),
  79 temp1:ev(eqlam,sub,diff),
  80 print("do you know apriori that some taylor coefficients are 0"),
  81 nullans:read("y,n")
  82 )$
  83
  84
  85 g_poly(i,j):=block(
  86 /*this is a function to determine a particular taylor coefficient of the
  87 bifurcation equation g(amp,lam) = 0. it requires kowledge about the taylor
  88 coefficients of w(amp,lam). this knowledge is stored in the list database*/
  89        ls_list:cons([k=i,l=j],ls_list),
  90        if nullans = y then (
  91                             zeroans:read("is g_poly(",i,",",j,")identically
  92                                                 zero, y/n"),
  93                             if zeroans = y then return(bifeq[i,j]:0)),
  94        temp2:diff(temp1,amp,i,lam,j),
  95        derivsubst:true,
  96        /*set the derivatives diff(w,amp,i,lam,j) to zero*/
  97        temp3:subst(difnull(i,j),temp2),
  98        /*enter all information in database*/
  99        d_base_length:length(database),
 100        for ii thru d_base_length do
 101                 temp3:ev(subst(database[d_base_length+1-ii],temp3),diff),
 102        derivsubst:false,
 103        temp4:expand(ev(temp3,amp=0,lam=0,wnull,integrate)),
 104        /*project onto afun*/
 105        temp5:ratsimp(temp4.afun),
 106        bifeq[i,j]:integrate(trigreduce(temp5),xvar,0,len))$
 107
 108
 109
 110
 111 w_poly(i,j):=block(
 112 /* this function allows to determine iteratively any particular taylor
 113 coefficient of the function w(amp,lam). it returns a differential equation
 114 for the particular coefficient of interest (called zw1,zw2...). this d.e.
 115 is solved via solve_ode and the result is fed into database from feedw*/
 116        if nullans = y then (
 117                     zeroans:read("is diff(w(amp,",i,",lam,",j,") identically
 118                                         zero, y/n"),
 119                     if zeroans = y then (
 120                                  addbase:difnull(i,j),
 121                                  database:append(database,addbase),
 122                                  return(addbase))),
 123        wmax_diff:map(lambda([u],'diff(u,amp,i,lam,j) = concat(z,u)),w),
 124        temp2:diff(temp1,amp,i,lam,j),
 125        derivsubst:true,
 126        /*rename the derivatives diff(w,amp,i,lam,j) to zw */
 127        temp3:subst(wmax_diff,temp2),
 128        /*enter all information in stored in database*/
 129        d_base_length:length(database),
 130        for ii thru d_base_length do
 131                 temp3:ev(subst(database[d_base_length+1-ii],temp3),diff),
 132        derivsubst:false,
 133        temp4:ev(temp3,amp=0,lam=0,wnull,integrate),
 134        /*this is the projection q onto the range of the
 135          linear differential operator in the problem*/
 136        temp5:integrate(ev(temp4,zwnull).afun,xvar,0,len),
 137        temp6:temp4-temp5/norm*cfun,
 138        temp7:trigreduce(temp6),
 139        /*the set of o.d.e.'s to solve */
 140        w_de:expand(temp7),
 141        temp8:ev(w_de,vfun(zwlist,0)),
 142        /*if the particular solution of w_de is zero then w=0 */
 143        if nullist=temp8 then ( addbase:difnull(i,j),
 144                                database:append(database,addbase),
 145                                return(addbase)),
 146        temp9:solve_ode(temp8),
 147        feedw(temp9)
 148        )$
 149
 150
 151 solve_ode(exp):=(
 152 /*this function solve the d.e. w_de by a fourier mode ansatz*/
 153        trigfun:[],
 154        const:false,
 155        for i thru n do (
 156        /*determine the fourier modes*/
 157                     trig1:exp[i],
 158                     if trig1 # 0 then  (
 159                                  trig2:apply1(trig1,sinnull,cosnull),
 160                                  if trig2 # 0 then (
 161                                               const:true,
 162                                               trig1:trig1-trig2),
 163                     trigfun:append(find_trig(trig1),trigfun))),
 164        sol1:delete(dummy,setify(trigfun)),
 165        /*make an ansatz*/
 166        ansatz:makelist(sum(am[i,j]*sol1[i],i,1,length(sol1)),j,1,n),
 167        sol2:ev(w_de,map("=",zwlist,ansatz),diff),
 168        sol3:makelist(ratcoef(sol2,i),i,sol1),
 169        eqlist:[],
 170        for i thru length(sol3) do eqlist:append(eqlist,sol3[i]),
 171        varlist:[],
 172        for i thru n do
 173                 for j thru length(sol1) do varlist:cons(am[j,i],varlist),
 174        /*find the amplitudes of the fourier modes*/
 175        sol4:solve(eqlist,varlist),
 176        /*solve for the constant fourier mode if necessary*/
 177        cansatz:0,
 178        if const = true then (cansatz:makelist(concat(con,i),i,1,n),
 179                              csol1:ev(w_de,map("=",zwlist,cansatz),diff),
 180                              csol2:apply1(csol1,sinnull,cosnull),
 181                              csol3:solve(csol2,cansatz)),
 182        ev(ansatz+cansatz,sol4,csol3))$
 183
 184
 185 feedw(exp):=(
 186 /* this function allows the result of the ode-solver to be entered into
 187    the list database. it checks for orthogonality to the critical adjoint
 188    eigenfunction and removes  solutions of the homogeneous equation (i.e.
 189    nonorthogonal terms)*/
 190        f2:integrate(exp.afun,xvar,0,len),
 191        if ratsimp(f2)=0
 192                 then
 193                      addbase:map("=",makelist('diff(w[u],amp,i,lam,j),u,1,n),
 194                                                         exp)
 195                 else (ortho_result:ratsimp(exp- f2/norm*cfun),
 196                      addbase:map("=",makelist('diff(w[u],amp,i,lam,j),u,1,n),
 197                                                 ortho_result)),
 198        database:append(database,addbase),
 199        print(addbase))$
 200 /*collect all information stored in bifeq and assemble the bifurcation
 201 equation*/
 202 g_eq():=  sum(ev(amp^k*lam^l/k!*bifeq[k,l],ls_list[n]),n,1,length(ls_list))$
 203
 204
 205 setify(list):=(/*transforms a list into a set, i.e. removes duplicates*/
 206         set:[list[1]],
 207         for i :2 thru length(list) do
 208                if not member(list[i],set) then
 209                         set:cons(list[i],set),
 210         set)$
 211
 212 find_trig(exp):=(/*finds the fourier modes*/
 213         f_a1:args(exp+dummy),
 214         f_a2:apply1(f_a1,sinfind,cosfind)
 215         )$
 216
 217 /*auxiliary functions */
 218 matchdeclare([dummy1,dummy2],true)$
 219 defrule(cosfind,dummy1*cos(dummy2),cos(dummy2))$
 220 defrule(sinfind,dummy1*sin(dummy2),sin(dummy2))$
 221 defrule(cosnull,cos(dummy1),0)$
 222 defrule(sinnull,sin(dummy1),0)$
 223
 224
 225 vfun(list,value):=map(lambda([u],u=value),list)$
 226 difnull(i,j):=map(lambda([u],'diff(u,amp,i,lam,j)=0),w) $