ibp.mpl Version 0.2

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. An example run would look like 


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! 

#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;
	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.2.";

	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
							newexpr := newexpr - op(i,expr2) - diff(thecoeff,dparam)*diff(function,dparam $ power-1);
						elif power > 0 then
							newexpr := newexpr - op(i,expr2) - diff(thecoeff,dparam)*function;
						else
							print("Problem! You must lower terms with a positive integer power!");
							RETURN(expr2);
						fi:
					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
						newexpr :=  - diff(thecoeff,dparam)*diff(function,dparam $ power-1);
					elif power > 0 then
						newexpr :=  - diff(thecoeff,dparam)*function;
					else
						print("Problem! You must lower terms with a positive integer power!");
						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(1,op(j,term));	
							#counter := 0;
							for counter from 2 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. If you are aware of this, pay no heed to me.");
									else
										thecoeff := coeff(term, op(j,term));
										newexpr := newexpr - op(i,expr2) - diff(thecoeff,dparam)*op(1,op(j,term));
									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(1,op(j,term));	
						#counter := 0;
						for counter from 2 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. If you are aware of this, pay no heed to me.");
								else
									thecoeff := coeff(term, op(j,term));
									newexpr := - diff(thecoeff,dparam)*op(1,op(j,term));
								fi:
							fi:
						end do:
					fi:
				fi:
			end do:
		fi:
	fi:

	RETURN(simplify(newexpr,size));
end: