Make sure MATLAB does not recalculate symbolic exp

Question

I am building (my first...) MatLab program, it needs to differentiate an equations symbolically and then use this solution many many times (with different numeric inputs).

I do not want it to recalculate the symbolic differentiation every time it needs to put in a new set of numeric values. This would probably greatly add to the time taken to run this program (which - given its nature, a numeric optimiser, will probably already be hours).

My question is how can I structure my program such that it will not recalculate the symbolic differentiation?

The class in question is:

function [ result ] = GradOmega(numX, numY, numZ, numMu)
syms x y z mu
omega = 0.5*(x^2+y^2+z^2) + (1-mu)/((x+mu)^2+y^2+z^2)^0.5 + mu/((x+mu-1)^2+y^2+z^2)^0.5;
symGradient = gradient(omega);
%//Substitute the given numeric values back into the funtion
result = subs(symGradient, {x,y,z,mu}, {numX, numY, numZ, numMu});
end

I know that I could just symbolically calculate the derivative and then copy-paste it into the code e.g.

gradX = x + ((2*mu + 2*x)*(mu - 1))/(2*((mu + x)^2 + y^2 + z^2)^(3/2)) - (mu*(2*mu + 2*x - 2))/(2*((mu + x - 1)^2 + y^2 + z^2)^(3/2));
gradY = y - (mu*y)/((mu + x - 1)^2 + y^2 + z^2)^(3/2) + (y*(mu - 1))/((mu + x)^2 + y^2 + z^2)^(3/2);
gradZ = z - (mu*z)/((mu + x - 1)^2 + y^2 + z^2)^(3/2) + (z*(mu - 1))/((mu + x)^2 + y^2 + z^2)^(3/2);

But then my code is a bit cryptic, which is a problem in a shared project. There is a related query here: http://uk.mathworks.com/matlabcentral/answers/53542-oop-how-to-avoid-recalculation-on-dependent-properties-i-hope-a-mathwork-developer-could-give-me-a But I'm afraid I couldn't follow the code. Also I am much more familiar with Java and Python, if that helps explain anything.

Answer 1

You could wrap your function into some kind of Function-Factory, which does not return numerical results, but a function that can be evaluated:

(I had to replace the call syms with sym('mu'), because for some reason it kept calling a mutools function in line omega = .... I did also change the call to gradient to make sure the arguments are in correct order, and mu will be treated as constant.)

function GradOmega = GradOmegaFactory()
x = sym('x');
y = sym('y');
z = sym('z');
mu = sym('mu');
omega = 0.5*(x^2+y^2+z^2) + (1-mu)/((x+mu)^2+y^2+z^2)^0.5 + mu/((x+mu-1)^2+y^2+z^2)^0.5;
symGradient = gradient(omega,{'x','y','z'});
GradOmega = matlabFunction(symGradient, 'vars', {'x','y','z','mu'});
end

Then you would call it via:

GradOmega = GradOmegaFactory();
result1 = GradOmega(numX1, numY1, numZ1, numMu1);
result2 = GradOmega(numX2, numY2, numZ2, numMu2);
result3 = GradOmega(numX3, numY3, numZ3, numMu3);
...

Even better:

You could go even more fancy and use a wrapper function GradOmega which builds such a function inside and makes it persistent, to get the same interface you had with your initial approach. The first time you call the function GradOmega the symbolic expression is evaluated, but on each consecutive call you will only have to evaluate the generated function handle, which means it should be nearly as fast as if you hard-coded it.

function result = GradOmega(numX, numY, numZ, numMu)
persistent numericalGradOmega;
if isempty(numericalGradOmega)
    numericalGradOmega = GradOmegaFactory();
end
result = numericalGradOmega(numX, numY, numZ, numMu);
end

Use this like you would use your original version

result = GradOmega(numX, numY, numZ, numMu);

Just copy and paste both functions into a single GradOmega.m file. (GradOmega should be the first function in the file.)

---

Another tip: You can even evaluate this function using vectors. Instead of calling GradOmega(1,2,3,4) and GradOmega(5,6,7,8) afterwards, you can save the time overhead via the call GradOmega([1,5], [2,6], [3,7], [4,8]) using row vectors.

Yet another tip: To clean up your code even more, you could also put the first lines into a separate symOmega.m file.

function omega = symOmega()
x = sym('x');
y = sym('y');
z = sym('z');
mu = sym('mu');
omega = 0.5*(x^2+y^2+z^2) + (1-mu)/((x+mu)^2+y^2+z^2)^0.5 + mu/((x+mu-1)^2+y^2+z^2)^0.5;

This way you don't have to have a copy of this symbolic expression in every file you use it. This can be beneficial if you also want to evaluate Omega itself, as you then can make use of the same Factory-approach listed in this answer. You would end up with the following files: symOmega.m, Omega.m and GradOmega.m, where only the file symOmega.m has the actual mathematical formula and the other two files make use of symOmega.m.