matlabFunction() generating intermediate terms that place the integrand outside integral()

I'm using the symbolic math toolbox to calculate a couple very large matrices containing integrals and then converting them to regular (non symbolic) functions via matlabFunction(). The problem is that MATLAB generates a bunch of intermediate terms (~200), most of which contain the variable I'm integrating over, which breaks the script. For exmaple:

M = integral(@(s)f(et1,et2,...,s)) ,<-- et1,et2,... = g(s)

I'm not sure if that's the absolute best way to put it, but I somehow need to find a way to either suppress all the extra terms or roll the extra terms containing s back into the bigger equation, since s isn't defined outside the integral.

So far I've tried:

• Playing with the settings for matlabFunction() - turning off optimization makes it generate even more extra terms, which doesn't help
• Importing the data in Excel and trying to manually substitute it in, which I don't think is a great idea because it involes about an hour of setup and would need to be redone every time I change the code at all (although that's probably because I'm not good enough at VBA to make it into a macro)
• Creating a global symbolic variable s and using that, which makes the calculation work but also kills my speed, making the 1/5 of the timespan I tested take over 6 hours to calculate

Does anyone have any ideas how to fix this? Is there a setting in MATLAB I'm missing?

[Paul]

Hi Alex,

Can you show a simple example with actual code that illustrates the problem, or at least illustrates the actual workflow? I'm not sure how to interpret this:

M = integral(@(s)f(et1,et2,...,s)) ,<-- et1,et2,... = g(s)

[Torsten]

@Alex B comment moved here:

syms q0 tht s D d L m phi g t
nSimplify = 1;
q = [q0; tht];
alpha = s*tht + (s^2 * q0)/2;
x = D*d*cos(alpha) - L*int(sin(alpha),s,0,s);
y = D*d*sin(alpha) + L*int(cos(alpha),s,0,s);
dq = real([gradient(x,q),gradient(y,q)]);
fprintf(" building inertia matrix...\n")
M = int(int(m*dq'*dq,d,-1/2,1/2),s,0,1)
matlabFunction(M,"file",inertiaMat)

This is an example of what I mean (it takes ~a minute to run for me); M will end up containing integrals in d and s with terms containing d and s inside them, but when matlabFunction returns the the file at the end it will generate intermediate terms containing s and d, find those, realize s and d aren't defined (because I'm outside th integral call) and crash.

Answer 1

[Paul]

A couple of suggestions on the code.

As a general rule IMO, it's best to specify all relevant assumptions on the variables. For example, I think that all of these variables are real.

syms q0 tht s D d L m phi g t real

And if we know that m (mass?) is positive:

assumeAlso(m,'positive')

Use assumeAlso to reflect any other assumptions on the variables.

nSimplify = 1;
q = [q0; tht];
alpha = s*tht + (s^2 * q0)/2;

I'm always surprised that Matlab allows the limits of integration to be expressed in terms of the variable of integration. I'll assume that this code is actually doing what's desired, but it might be better to use subs() on alpha to integrate wrt a dummy variable of integration

x = D*d*cos(alpha) - L*int(sin(alpha),s,0,s)
y = D*d*sin(alpha) + L*int(cos(alpha),s,0,s)

Note sue why orginal code had a need to use real() on the next line.

dq = [gradient(x,q),gradient(y,q)];
% fprintf(" building inertia matrix...\n")

Use .' (note the dot). Using just ' means complex conjugate transpose, which isn't what's wanted here I don't believe.

Also, it's good practice to use the parentheses to ensure symmetry.

M = int(int(m*(dq.'*dq),d,-1/2,1/2),s,0,1)

@Walter Roberson stresses to use Optimize/false in creating the file, so do that here.

matlabFunction(M,"file",'inertiaMat','Optimize',false)

Having said all of that, I see exactly the problem, which is that inertiaMat.m contains a line like

et2 = (tht.*sqrt(pi).*fresnels(1.0./sqrt(pi).*(tht./2.0+(q0.*s)./2.0).*sqrt(1.0./q0).*2.0).*sin(tht.^2./(q0.*2.0)).*sqrt(1.0./q0))./q0;

which in this context should be:

et2 = @(s) (tht.*sqrt(pi).*fresnels(1.0./sqrt(pi).*(tht./2.0+(q0.*s)./2.0).*sqrt(1.0./q0).*2.0).*sin(tht.^2./(q0.*2.0)).*sqrt(1.0./q0))./q0;

and then all uses of et2 in integral() should be converted to et2(s). Same for all of the other et variables.

I don't know why matlabFunction does what it does; maybe there's a sympref that controls this behavior. Can't think of a work around other than editing inertiaMat.m after it's created.

[Alex B]

Thank you so much for all the good info! I'll look through all my symprefs and see if there's somethign there

(I'll accept this answer depending if Walter gets back about this being a bug -- moving up a version to fix this is an option)

[Alex B]

Sorry for the super late necro-reply but I can confirm that MATLAB r2022a doesn't seem to exhibit this bug and generates integrals properly

Answer 2

[Walter Roberson]

You are correct, that is a known bug in your release, and there is no solution in your release.

[Paul]

I'm using 2021b and still seeing it. Was it fixed in 2022a to your knowledge?

[Alex B]

I'm not seeing anything about it in the 21b/22a release notes so I'm going to assume no here