## ibp.mpl Version 0.5

Here is a utility (ibp) that allows you to integrate an expression (not an integral) by parts in a manner that could cause permanent damage to any nearby mathematicians. It neglects boundary terms and returns an expression.

The primary utility is the removal of high-order derivatives from an expression in, for example, the variation of an action. If you wish to show that an action leads to second-order equations of motion, sometimes you must perform integration by parts before a functional variation in order to see that the field equations are indeed second-order.

It will also attempt to integrate a function whose integration by parts procedure is recursive. That is, if it detects that function --> k*function + terms when integrated by parts (where k is some constant), it will attempt to solve the resulting equation.

In addition, it features an option, tryhard, which is a boolean true/false that will tell it to attempt to force integration by parts on expressions with non-unity powers of derivatives. In that case, even though it will not be possible to be rid of the derivative, the procedure will still integrate it by parts.

### Example


expr := -2*g(r)^2*f(r)*diff(N(r),r$2) - g(r)^4*N(r)*diff(f(r),r$2)
+ f(r)*g(r)*diff(N(r),r$3); newexpr := ibp(expr, 3, r); newexpr; (-g(r)*diff(f(r),r)-2*f(r)*(g(r)^2+1/2*diff(g(r),r)))* diff(diff(N(r),r),r)-g(r)^4*N(r)*diff(diff(f(r),r),r)   No third derivatives! Hooray! ### Source Code  #This maple procedure will integrate an algebraic function by parts, neglecting #boundary terms #First, the helper procedures _ibp_solveeqn := proc(outexpr,oldterm,ibpcoeff,derivterm,dparam) local newterm, counter, myinteger; #outexpr is the full expression that is returned by this procedure. #oldterm is the original term that we are integrating by parts #ibpcoeff is the coefficient which we treat as V in VdU #derivterm is the term to integrate out...it needs to be in the form #diff(f(x),x$n)
#where n is some integer!
newterm := expand(diff(ibpcoeff,dparam)*op(1,derivterm));
if is(op(0,newterm) = +) then
for counter from 1 to nops(newterm) do:
if is(op(counter,newterm)/oldterm,real) then
#If there exists a term in the IBP expansion equal to the original term up to a constant,
#solve the equation to get a closed-form solution for oldterm.
myinteger := op(counter,newterm)/oldterm;
newterm := ( newterm - myinteger*oldterm)/(1+myinteger);
break;
fi:
end do:
fi:
RETURN(outexpr - oldterm - newterm);
end proc:

_ibp_specifiedfn := proc(functions,newexpr,term,dparam,power,tryhard)
local function, j, localexpr,exitparam, tmp, thecoeff, derivterm;
exitparam := false;
localexpr:=newexpr;
for function in functions do
if (has(term,diff(function,dparam $power)) and not(has(term,diff(function,dparam$ power+1))) ) then
if is(op(0,term)=*) then
for j from 1 to nops(term) do
tmp := op(j,term);
#What if tmp = diff? i.e. expr := 2*diff(f(r),r) + diff(f(r),r)
if (has(tmp,diff(function,dparam $power)) and not(has(tmp,diff(function,dparam$ power+1))) ) then
if is(op(0,tmp) = ^) then
thecoeff := coeff(term,diff(function,dparam $power)^op(2,tmp))*diff(function,dparam$ power)^(op(2,tmp)-1);
if is(tryhard=false) then
exitparam := true;
fi:
break;
elif is(op(0,tmp) = diff) then
thecoeff := coeff(term,diff(function,dparam \$ power));
break;
else
print("Error - could not determine power of term.");
RETURN(localexpr);
fi:
fi:
end do:
if is(exitparam = true) then
break;
else
derivterm :=term/thecoeff;
localexpr := _ibp_solveeqn(localexpr,term,thecoeff,derivterm,dparam);
break;
fi:
else
#Do nothing - this is a lone term!
#Unless it is a power of derivatives!
#In this case we can get a higher power derivative, (if we tryhard)
if(op(0,term)=^) then
tmp := op(1,term);
#What if tmp = diff? i.e. expr := 2*diff(f(r),r) + diff(f(r),r)
if is(op(0,tmp) = diff) then
if ((tryhard)) then
thecoeff := tmp^(op(2,term)-1);
derivterm :=term/thecoeff;
localexpr := _ibp_solveeqn(localexpr,term,thecoeff,derivterm,dparam);
else
RETURN(localexpr);
fi:
else
print("Error - could not parse the term",term,".");
RETURN(localexpr);
fi:
fi:
fi:
fi:
end do:
RETURN(localexpr);
end proc:

ibp := proc(expr,power,dparam,functions:=[],tryhard:=false)
local expr2,newexpr,exprlength,termlength,term,thecoeff,i,j,tmp,counter,derivterm;
description "This utility integrates an expression by parts in the manner most likely to make a mathematician cry, picking out the terms of power \'power\' and differentiating with respect to the other terms. Boundary terms are always neglected. Arguments are: the expression, the power of derivative you wish to lower, the variable with respect to which you differentiate, a list of the functions which possess the higher derivatives which you wish to integrate (optional), and a boolean determining whether we integrate expressions inside of powers or not (optional - default no). If you do not specify the final option, I will assume you want to integrate all function that have a derivative of order \'power\'. Written by W. Brenna, Aug 14 2012. Version 0.5 - modified January 27, 2013.";

expr2 := expand(expr);
exprlength := nops(expr2);

newexpr := expr2;

#Make sure we have a reasonable power.
if not(power > 0 and is(power,integer)) then
print("Error! Currently fractional or negative powers are unsupported.");
RETURN(expr2);
fi:

if not(is(functions=[])) then
if is(op(0,expr2)=+) then
for i from 1 to exprlength do
term := op(i,expr2);
newexpr := _ibp_specifiedfn(functions,newexpr,term,dparam,power,tryhard):
end do:
else
term := expr2;
newexpr := _ibp_specifiedfn(functions,newexpr,term,dparam,power,tryhard):
fi:
else
#populate the list of functions to be "all functions", if none were specified
if is(op(0,expr2) = +) then
for i from 1 to exprlength do
term := op(i,expr2);
termlength := nops(term);
for j from 1 to termlength do
if is(op(0,op(j,term)) = diff) then
if is(op(2,op(j,term)) = dparam) then
tmp := op(j,term);
for counter from 1 to power+1 while is(op(0,tmp) = diff) do
tmp := op(1,tmp);
if is(counter = power) then
if is(op(0,tmp) = diff) then
print("Warning: you have a derivative higher than the power you specified. It was ignored. If you are aware of this, pay no heed to me.");
else
#bonus code to perform infinite integration by parts!
thecoeff := coeff(term, op(j,term));
derivterm := term/thecoeff;
newexpr := _ibp_solveeqn(newexpr,term,thecoeff,derivterm,dparam);
fi:
fi:
end do:
fi:
fi:
end do:
end do:
elif is(op(0,expr2) = ^) then
#This is a lone term.
print("This is a lone term. You need to specify tryhard=true and explicitly give functions in order to solve this.");

else
term := expr2;
termlength := nops(term);
for j from 1 to termlength do
if is(op(0,op(j,term)) = diff) then
if is(op(2,op(j,term)) = dparam) then
tmp := op(j,term);
for counter from 1 to power+1 while is(op(0,tmp) = diff) do
tmp := op(1,tmp);
if is(counter = power) then
if is(op(0,tmp) = diff) then
print("Warning: you have a derivative higher than the power you specified. It was ignored. If you are aware of this, pay no heed to me.");
else
thecoeff := coeff(term, op(j,term));
derivterm := term/thecoeff;
newexpr := _ibp_solveeqn(newexpr,term,thecoeff,derivterm,dparam);
fi:
fi:
end do:
fi:
fi:
end do:
fi:
fi:
RETURN(simplify(newexpr,size));
end: