perturb := proc () # #Expands a set/list of equations in terms of variables in #in a variable list to specified order. Perturbation takes #form of: # x -> x[0] + epsilon*x[1] # where # x[1](r,t) = x[1](r)*exp(sigma*t) # #Takes as input # eqslist --- set or list of equations to perturb # varlist --- list of variable to perturb from which # it builds a list of expansions to subs. # optional: # parameter --- expansion parameter; Default: epsilon # order --- order to which the perturbation should # be carried; Default: 1 # Sample Usage: # ( from Holmes, Introduction to Perturbation Methods, pg. 275 ) # # read `pr.ms`; # testvars:= {mu}; # atest:=mu_rr + 2*beta*mu_r + F(mu) = 0: # test := mu_tt + mu_t = mu_rr + lambda*mu + f(mu): # a1:={test,atest}: # perturb(a1,testvars,delta,1); # # local result1, result2, result3, result4, eqn, order, parameter, eqslist, expansionlist,varlist, var, elem, elem2, elem3, elem4, elem5, elem6, elem7, elem8; if nargs<4 then order := 1 else order := args[4] fi; if nargs<3 then parameter := epsilon else parameter := args[3] fi; if nargs<2 then ERROR(`Need at least `) fi; if type(args[1],equation) then eqslist := {args[1]}; elif type(args[1],set) or type(args[1],list) then eqslist := args[1]; else ERROR(`You must input an equation, or a list or set of equations`) fi; if type(args[2],string) then varlist := {args[2]}; elif type(args[2],set) or type(args[2],list) then varlist := args[2]; else ERROR(`You must input a variable, or a list or set of variables`) fi; # build expansionlist from variable list expansionlist := []; for var in varlist do elem := var = ``.var.``[0] + parameter*``.var.``[1]; elem2 := ``.var._r = ``.var._r.``[0] + parameter*``.var._r.``[1]; elem3 := ``.var._rr = ``.var._rr.``[0] + parameter*``.var._rr.``[1]; elem4 := ``.var._t = parameter*sigma*``.var.``[1]; elem5 := ``.var._tt = parameter*sigma**2*``.var.``[1]; elem6 := ``.var._ttr = parameter*sigma**2*``.var._r.``[1]; elem7 := ``.var._tr = parameter*sigma*``.var._r.``[1]; elem8 := ``.var._rt = parameter*sigma*``.var._r.``[1]; expansionlist := [op(expansionlist), elem,elem2,elem3,elem4,elem5,elem6,elem7, elem8]; od; # done building "expansionlist" result1:= simplify( subs( expansionlist, eqslist ) ); result4:= []; for eqn in result1 do result2:= simplify(convert(series(lhs(eqn),parameter,order+1),polynom) = convert(series(rhs(eqn),parameter,order+1),polynom) ); result3:= simplify( subs(parameter=0,result2) ); result4:= [op(result4),simplify( (result2-result3)/parameter)]; od; result4 end: takederiv := proc(eqn) local su1, su2, su3, eqn2, eqn3, eqn4, eqn5; su1:={ omega[0] = omega[0](r), lambda[0] = lambda[0](r), nu[0] = nu[0](r), Gamma[0] = Gamma[0](r), psi[0] = psi[0](r), omega[1] = omega[1](r), lambda[1] = lambda[1](r), nu[1] = nu[1](r), Gamma[1] = Gamma[1](r), psi[1] = psi[1](r), omega = omega(r), lambda = lambda(r), nu = nu(r), Gamma = Gamma(r), psi = psi(r), omega_r[1] = omega_r[1](r), lambda_r[1] = lambda_r[1](r), nu_r[1] = nu_r[1](r), Gamma_r[1] = Gamma_r[1](r), psi_r[1] = psi_r[1](r), omega_r[0] = omega_r[0](r), lambda_r[0] = lambda_r[0](r), nu_r[0] = nu_r[0](r), Gamma_r[0] = Gamma_r[0](r), psi_r[0] = psi_r[0](r), omega_r = omega_r(r), lambda_r = lambda_r(r), nu_r = nu_r(r), Gamma_r = Gamma_r(r), psi_r = psi_r(r), omega_rr[1] = omega_rr[1](r), lambda_rr[1]= lambda_rr[1](r), nu_rr[1] = nu_rr[1](r), Gamma_rr[1] = Gamma_rr[1](r), psi_rr[1] = psi_rr[1](r), omega_rr[0] = omega_rr[0](r), lambda_rr[0]= lambda_rr[0](r), nu_rr[0] = nu_rr[0](r), Gamma_rr[0] = Gamma_rr[0](r), psi_rr[0] = psi_rr[0](r), omega_rr = omega_rr(r), lambda_rr = lambda_rr(r), nu_rr = nu_rr(r), Gamma_rr = Gamma_rr(r), psi_rr = psi_rr(r) }: eqn2 := subs( su1, eqn ); eqn3 := simplify(diff(eqn2,r)); su2:={ diff(omega[0](r),r) = omega_r[0], diff(lambda[0](r),r) = lambda_r[0], diff(nu[0](r),r) = nu_r[0], diff(Gamma[0](r),r) = Gamma_r[0], diff(psi[0](r),r) = psi_r[0], diff(omega[1](r),r) = omega_r[1], diff(lambda[1](r),r) = lambda_r[1], diff(nu[1](r),r) = nu_r[1], diff(Gamma[1](r),r) = Gamma_r[1], diff(psi[1](r),r) = psi_r[1], diff(omega(r),r) = omega_r, diff(lambda(r),r) = lambda_r, diff(nu(r),r) = nu_r, diff(Gamma(r),r) = Gamma_r, diff(psi(r),r) = psi_r, diff(omega_r[0](r),r) = omega_rr[0], diff(nu_r[0](r),r) = nu_rr[0], diff(lambda_r[0](r),r) = lambda_rr[0], diff(Gamma_r[0](r),r) = Gamma_rr[0], diff(psi_r[0](r),r) = psi_rr[0], diff(omega_r[1](r),r) = omega_rr[1], diff(nu_r[1](r),r) = nu_rr[1], diff(lambda_r[1](r),r) = lambda_rr[1], diff(Gamma_r[1](r),r) = Gamma_rr[1], diff(psi_r[1](r),r) = psi_rr[1], diff(omega_r(r),r) = omega_rr, diff(nu_r(r),r) = nu_rr, diff(lambda_r(r),r) = lambda_rr, diff(Gamma_r(r),r) = Gamma_rr, diff(psi_r(r),r) = psi_rr }: eqn4 := simplify(subs(su2,eqn3)); su3:={ omega[0](r) = omega[0], lambda[0](r) = lambda[0], nu[0](r) = nu[0], Gamma[0](r) = Gamma[0], psi[0](r) = psi[0], omega[1](r) = omega[1], lambda[1](r) = lambda[1], nu[1](r) = nu[1], Gamma[1](r) = Gamma[1], psi[1](r) = psi[1], omega(r) = omega, lambda(r) = lambda, nu(r) = nu, Gamma(r) = Gamma, psi(r) = psi, omega_r[0](r) = omega_r[0], lambda_r[0](r) = lambda_r[0], nu_r[0](r) = nu_r[0], Gamma_r[0](r) = Gamma_r[0], psi_r[0](r) = psi_r[0], omega_r[1](r) = omega_r[1], lambda_r[1](r) = lambda_r[1], nu_r[1](r) = nu_r[1], Gamma_r[1](r) = Gamma_r[1], psi_r[1](r) = psi_r[1], omega_r(r) = omega_r, lambda_r(r) = lambda_r, nu_r(r) = nu_r, Gamma_r(r) = Gamma_r, psi_r(r) = psi_r }: eqn5 := simplify( subs( su3, eqn4) ); end: reverse_subs := proc(eqlist) local temp1, temp2, rlist; rlist := {}; for temp1 in eqlist do temp2 := rhs(temp1) = lhs(temp1); rlist := rlist union {temp2}; od; rlist; end: latex_eq := proc( eqn, doit, file, add_or_overwrite) if (doit) then if (add_or_overwrite = `add`) then appendto(file); elif (add_or_overwrite = `over`) then writeto(file); else print(`not recognized option. Defaulting to add`); appendto(file); fi; lprint(`$$`); latex(eqn); lprint(`$$`); writeto(terminal); fi: end: simplecollect := proc(expr, varlist) # # simplecollect( , # ) # # Purpose: spits out list of coefficients of the terms in the # variable list (in the order input) in the same # fashion as collect, but also simplifies these # terms independently of each other. # # Returns: a list of coefficients in the same order as the input # variable list, *but* of length nops(varlist)+1 where # the last term contains all terms lacking a term in # the variable list. # # Note: if any cross-terms appear in the input, they do # not appear in the output. Hence, this should be used # only for an independent list of variables. # # Sample usage: # # simplecollect(a*x+10*y+x+exp(nu)*z+20 + x*y,[x,y,z]); # returns: # [a+1,10,exp(nu),20] # # local A, i, j, tsub, sub, all,therest; A := []; all := {}; if ( nops(varlist) = 0 ) then A := expr; else for j from 1 to nops(varlist) do tsub := {}; for i from 1 to nops(varlist) do sub := { op(i,varlist) = 0 }; if ( i <> j) then tsub := tsub union sub; fi; if ( j = 1) then all := all union sub; fi; od; therest:= simplify(subs(all,expr)); A := [op(A), simplify( (subs( tsub, expr )-therest ) /op(j,varlist) ) ]; od; A := [op(A), therest]; fi; end: