ibp.mpl Version 0.3
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.
This edition contains a new feature - it will 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.
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
ibp := proc(expr,power,dparam,functions:=[])
local expr2,newexpr,exprlength,termlength,function,term,thecoeff,i,j,tmp,counter,counter2,myinteger,newterm;
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, and the variable with respect to which you differentiate, and a list of the functions which possess the higher derivatives which you wish to integrate (optional). 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.3 - modified Aug 21, 2012.";
expr2 := expand(expr);
exprlength := nops(expr2);
newexpr := expr2;
if not(is(functions=[])) then
if is(op(0,expr2)=`+`) then
for i from 1 to exprlength do
term := op(i,expr2);
for function in functions do
if (has(term,diff(function,dparam $ power)) and not(has(term,diff(function,dparam $ power+1))) ) then
thecoeff := coeff(term,diff(function,dparam $ power));
if power > 1 then
newterm := expand(diff(thecoeff,dparam)*diff(function,dparam $ power-1));
if is(op(0,newterm) = `+`) then
for counter2 from 1 to nops(newterm) do:
if is(op(counter2,newterm)/term,real) then
myinteger := op(counter2,newterm)/term;
newterm := ( newterm - myinteger*term)/(1+myinteger);
#print(newterm,myinteger*op(i,expr2),myinteger);
break;
fi:
end do:
fi:
#print(op(counter2,newterm),op(i,expr2));
newexpr := newexpr - op(i,expr2) - newterm;
elif power = 1 then
newterm := expand(diff(thecoeff,dparam)*function);
if is(op(0,newterm) = `+`) then
for counter2 from 1 to nops(newterm) do:
if is(op(counter2,newterm)/term,real) then
myinteger := op(counter2,newterm)/term;
newterm := ( newterm - myinteger*term)/(1+myinteger);
#print(newterm,myinteger*op(i,expr2),myinteger);
break;
fi:
end do:
fi:
#print(op(counter2,newterm),op(i,expr2));
newexpr := newexpr - op(i,expr2) - newterm;
else
print("Problem! Currently fractional or negative powers do not work with custom function. Try without specifying the function.");
RETURN(expr2);
fi:
#want a break function here to prevent continuing the cycle through functions on the same term!
fi:
end do:
end do:
else
term := expr2;
for function in functions do
if (has(term,diff(function,dparam $ power)) and not(has(term,diff(function,dparam $ power+1))) ) then
thecoeff := coeff(term,diff(function,dparam $ power));
if power > 1 then
newterm := expand(diff(thecoeff,dparam)*diff(function,dparam $ power-1));
if is(op(0,newterm) = `+`) then
for counter2 from 1 to nops(newterm) do:
if is(op(counter2,newterm)/term,real) then
myinteger := op(counter2,newterm)/term;
newterm := ( newterm - myinteger*term)/(1+myinteger);
#print(newterm,myinteger*op(i,expr2),myinteger);
break;
fi:
end do:
fi:
#print(op(counter2,newterm),op(i,expr2));
newexpr := -newterm;
elif power = 1 then
newterm := expand(diff(thecoeff,dparam)*function);
if is(op(0,newterm) = `+`) then
for counter2 from 1 to nops(newterm) do:
if is(op(counter2,newterm)/term,real) then
myinteger := op(counter2,newterm)/term;
newterm := ( newterm - myinteger*term)/(1+myinteger);
#print(newterm,myinteger*op(i,expr2),myinteger);
break;
fi:
end do:
fi:
#print(op(counter2,newterm),op(i,expr2));
newexpr := -newterm;
else
print("Problem! Currently fractional or negative powers do not work with custom function. Try without specifying the function.");
RETURN(expr2);
fi:
fi:
end do:
fi:
else
#populate the list of functions to be "all functions", if none were specified
printf("I'm assuming you want me to integrate all functions of order %g.\n",power);
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));
newterm := expand(diff(thecoeff,dparam)*op(1,op(j,term)));
if is(op(0,newterm) = `+`) then
for counter2 from 1 to nops(newterm) do:
if is(op(counter2,newterm)/term,real) then
myinteger := op(counter2,newterm)/term;
newterm := ( newterm - myinteger*term)/(1+myinteger);
#print(newterm,myinteger*op(i,expr2),myinteger);
break;
fi:
end do:
fi:
#print(op(counter2,newterm),op(i,expr2));
newexpr := newexpr - op(i,expr2) - newterm;
fi:
fi:
end do:
fi:
fi:
end do:
end do:
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));
newterm := expand(diff(thecoeff,dparam)*op(1,op(j,term)));
if is(op(0,newterm) = `+`) then
for counter2 from 1 to nops(newterm) do:
if is(op(counter2,newterm)/term,real) then
myinteger := op(counter2,newterm)/term;
newterm := ( newterm - myinteger*term)/(1+myinteger);
#print(newterm,myinteger*op(i,expr2),myinteger);
break;
fi:
end do:
fi:
#print(op(counter2,newterm),op(i,expr2));
newexpr := -newterm;
fi:
fi:
end do:
fi:
fi:
end do:
fi:
fi:
RETURN(simplify(newexpr,size));
end:
![wbrenna [\at/] uwaterloo [\dot/] ca](/wilson/includes/email.png)