• Home
• CV
• Research
• Computing
• Projects

physicsdiff.mpl Version 0.2

This procedure allows you to perform a first variation of some expression (presumed to be inside an integral, e.g. the Lagrangian density in an action). It will differentiate the expression treating the specified function as the differentiation parameter, with the additional feature that it automatically attempts to integrate first order derivatives of that function by parts.

To call the procedure, run: physicsdiff(expression,function,diffparam);

Example

expr := h(r)*g(r)*2 + h(r)*diff(h(r),r)^2*3 + diff(h(r),r):
physicsdiff(expr,h(r),r);
	2*g(r) + 3*(diff(h(r),r))^2 - 2*diff(3*h(r)*diff(h(r),r),r)

Source Code

#This maple procedure does a functional variation, much like diff does when used
#with(Physics):, but with the added feature that it can vary simple derivatives
#with respect to that function, ie. vary(c*diff(f(r),r)^2,f(r)) using
#integration by parts.

_physicsdiff2_term := proc(function,term,innerterm,dparam)
	local counter, thecoeff,tmp;

	if is(op(0,innerterm) = `^`) then
#In this case, the term is a power of a derivative of the function of interest
		thecoeff := coeff(term,diff(function,dparam)^(op(2,innerterm)));
		#print(term,thecoeff,innerterm,"isntdiff");
		RETURN(op(2,innerterm)*diff(thecoeff*diff(function,dparam)^(op(2,innerterm)-1),dparam));
	elif is(op(0,innerterm)=`diff`) then
#Here the function of interest is alone
		thecoeff := coeff(term,diff(function,dparam));
		#print(term,innerterm,"isdiff");
		#RETURN(diff(thecoeff,dparam)*function);
#We can't return the function with it! That is the part that gets differentiated away!
		RETURN(diff(thecoeff,dparam));
	else
		print("Error - could not determine power of term.");
		RETURN(0);
	fi:

end proc:


_physicsdiff_term := proc(expr,function,term,dparam)
	local counter,counter2,tmp,retterm;
	retterm := 0;

	if is(op(0,dparam)=`list`) then
		for counter from 1 to nops(dparam) do
			if (has(term,diff(function,dparam[counter])) and not(has(term,diff(function,dparam[counter] $ 2)))) then
				if is(op(0,term)=`*`) then
					for counter2 from 1 to nops(term) do
						tmp := op(counter2,term);
						if (has(tmp,diff(function,dparam[counter]))) then
							retterm := retterm - _physicsdiff2_term(function,term,tmp,dparam[counter]);
						fi:
					end do:
				elif is(op(0,term)=`^`) then
					tmp := term;
					if(has(tmp,diff(function,dparam[counter]))) then
						retterm := retterm - _physicsdiff2_term(function,term,tmp,dparam[counter]);
					fi:
				else
					print("Found a term that's just a lone derivative. Integration by parts failed.");
				fi:
				RETURN(retterm);
			elif has(term,diff(function,dparam[counter] $ 2)) then
				print("Error, found a derivative of order greater than 2. This term is killed.");
				RETURN(0);
			fi:
		end do:
		RETURN(0);
	else
		if (has(term,diff(function,dparam)) and not(has(term,diff(function,dparam $ 2)))) then
			if is(op(0,term)=`*`) then
				for counter from 1 to nops(term) do
					tmp := op(counter,term);
					if (has(tmp,diff(function,dparam))) then
						retterm := retterm - _physicsdiff2_term(function,term,tmp,dparam);
					fi:
				end do:
			elif is(op(0,term)=`^`) then
				print("Found a term that has powers.");
				tmp := term;
				if(has(tmp,diff(function,dparam))) then
					retterm := retterm - _physicsdiff2_term(function,term,tmp,dparam);
				fi:
			else
#There is nothing to be done - the term is just a lone derivative of "function" - it will be killed.
				print("Found a term that's just a lone derivative. Integration by parts failed.");
			fi:
			RETURN(retterm);
		elif has(term,diff(function,dparam $ 2)) then
			print("Error, found a derivative of order greater than 2. This term is killed.");
			RETURN(0);
		else
			RETURN(0);
		fi:
	fi:

end proc:




physicsdiff := proc(expr,function,dparam)
	local expr2,newexpr,exprlength,term,tmp,counter;
	description "This Maple procedure will perform a physicists' functional differentiation, neglecting boundary terms and integrating by parts to functionally vary an expression with respect to a function, even if it contains powers of first derivatives of that function. Written by W. Brenna, Sept 3 2012. Version 0.2.";

	if not(is(op(0,f(r)*g(r))=`*`)) then
		print("Error: Conflicting package detected! The `*` declaration has been redefined and results are therefore unreliable. Make sure you have loaded physicsdiff before running with(Physics).");
		RETURN(expr);
	fi:


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


	#newexpr := expr2;
	newexpr := 0;


	if is(op(0,expr2)=`+`) then
		for counter from 1 to exprlength do
			term := op(counter,expr2);
#This just does the first order integration by parts.
			newexpr := newexpr + _physicsdiff_term(newexpr,function,term,dparam):
		end do:
	else
		term := expr2;
		newexpr := newexpr + _physicsdiff_term(newexpr,function,term,dparam):
	fi:

#Finally, we do the regular differentiation by parts to build up the final
#terms.
	#RETURN(newexpr + frontend(diff,[expr2,function]));
#frontend is hopeless for this - non-integer powers make it cry.
	#print(Physics:-diff(expr2,function));
	RETURN(newexpr + Physics:-diff(expr2,function));

end proc:

---