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 <eqslist> <varlist>`) 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( <expression to be colected and simplified>,
# <ordered list of variables on which to collect> )
#
# 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: