matlabFunction() generating intermediate terms that place the integrand outside integral()
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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?
2 Commenti
Paul
il 21 Giu 2022
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
il 21 Giu 2022
@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.
Risposta accettata
Paul
il 21 Giu 2022
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.
2 Commenti
Più risposte (1)
Walter Roberson
il 21 Giu 2022
You are correct, that is a known bug in your release, and there is no solution in your release.
2 Commenti
Vedere anche
Categorie
Scopri di più su Assumptions in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!