My 6x6 symbolic Jacobian matrix is massive, and matlabFunction( ) is having trouble writing the corresponding numerical function file.

Hi there!

I currently have a 6x6 symoblic Jacobian using the Symbolic Math Toolbox. It is massive, when I take a peek at it using the Command Window. I then used matlabFunction( ) to convert this symbolic Jacobian to a numerical function. If I use matlabFunction's default optimize = true, the numerical function file is about 8,000 lines of code, and I get a message there that says the code is too complex to analyze, and that I should reduce the number of operations in my code. On the other hand, if I use matlabFunction's optimize = false option, the numerical function file is 59,506 lines of code, and I get a message there that says the file is too large. I then tried passing to matlabFunction( ), "Sparse", true, to try to get a sparse numerical Jacobian function file, but I still get the message in the numerical file that the code is too complex to analyze.

What can I do from here?

My goal is to successfully get a numerical Jacobian function file, then evaluate it to get a numerical Jacobian, and then find its eigenvalues using eig( ).

Thanks in advance,

[Torsten]

I currently have a 6x6 symoblic Jacobian using the Symbolic Math Toolbox.

If you already have the (6x6) Jacobian matrix in symbolic form, why don't you just use "subs" to substitute numerical values (your stationary points) for the y-vector and apply "eig" to the resulting numerical matrix ?

[Noob]

Hi Torsten!

I read something about garbage in garbage out. Are you saying I can simply bypass using matlabFunction and just evaluate my symbolic Jacobian numerically, and then find eigenvalues numerically? And use all numerical values for all physical parameters and forces laws etc? I had thought I couldn’t do this, and that I first had to write to a numerical file, and then do things numerically in a separate script file, which calls the numerical function file create by matlabFunction. So, those are my concerns. What do you think?

Thanks!

[Matt J]

The problem with evaluating symbolically subs() is that it will be very slow. If the calculation needs to be repeated, this could be a major bottleneck.

[Torsten]

I read something about garbage in garbage out.

If you have the Jacobian in symbolic form, you don't have garbage.

Since you only have a small number of points (or even a single point) where to evaluate the Jacobian, time is not a factor. So just proceed as I told you to proceed.

If J is your symbolic Jacobian, try

Jnum = double(subs(J,y,ynum))
eig(Jnum)

where ynum is the numerical point where you want to evaluate it and y is the name of your symbolic (6x1) solution vector.

That's all.

[Torsten]

If your ODE system is that complicated, you should completely avoid symbolic computations. Set up everything numerically right from the start - then you won't run into these problems you describe.

Maybe someone can help how to do this, but without at least knowing the mathematical model you try to simulate, it's almost impossible.

[Noob]

Hi Torsten!

Yes, working completely numerically is something I should consider in the near future perhaps.

What led me to work symbolically was because my equations are difficult in this sense:

Variables show up both on the LHS and RHS.

Symbolic math can handle this, such as x = x + 5 + y.

A numerical code couldn't handle this, and would say x is undefined.

The equal sign in a numerical code assigns numerical values to a variable.

If you want to share more of your thoughts, please feel free to do so.

Currently, I am separating the Jacobian and storing them into multiple variables, and I will then matlabFunction each of these smaller expressions.

I might also try using expand( ) followed by simplify( ), on all of the split-up Jacobian variables.

Thanks!

Answer 1

[Matt J]

Maybe just compute the Jacobian numerically, by finite differences? If the analytical expressions for the Jacobian are that massive, I don't see how you can be any more confident in the numerical accuracy of the analytical derivative than a finite difference derivative.

[Noob]

Hi Matt!

Is there any particular finite-difference method you’d recommend? Should I write it myself (following a book?) or use one that’s already developed, say, posted on the File Exchange?

Thanks!

[Matt J]

This is what I use,

https://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation?s_tid=ta_fx_results

[Noob]

Hi Matt!

After a few days of trying a lot of different things, I now think your advice is best to go with. I don't want to compromise my good workflow and make drastic changes to accommodate a massive symbolic Jacobian -- that seems to be the wrong direction to head in. So I will both look at the file you mentioned (it seems extremely popular) as well as try to learn to write a central-difference code myself.

Thanks again for your help!

Answer 2

[Walter Roberson]

Use children() to get a 6 x 6 cell array (or simply loop over the contents of the array); matlabFunction() each one seperately .

The time taken to use matlabFunction 'optimize' is roughly exponential in the size of the expression, so separating the expressions can help a lot.

If necessary, you can use children() on the individual expressions, if you are careful about how you put the pieces back together again.

[Noob]

Hi Walter!

Where do I write the children( ) command? After using the symbolic jacobian( ) command on my six first-order equations? Or do you mean write children( ) in the Command Window?

Thanks!

[Walter Roberson]

J = jacobian(EXPRESSION);
ch = children(J);
numch = numel(ch);
jh = cell(size(J));
for K = 1 : numch
    outfile = "jacob" + K + ".m";
    jh{K} = matlabFunction( ch{K}, 'file', outfile, 'optimize', true, 'vars', LIST_OF_VARIABLES);
end

Now jh will be a cell array of function handles that you can call with numeric values.

nJ = zeros(size(J));
for K = 1 : numch
    nJ(K) = jh{K}(var1, var2, var3, ...);
end

[Noob]

Hi Walter!

This is my idea inspired by you. Worth a try maybe? Here it is:

I think matlabFunction would be happy to manage these little blocks of symbolic code I just thought of:

1. In my symbolic code, use the symbolic differentiation command diff( ).

2. First equation, differentiate using diff with respect to each state variable. That means use diff six times. Now, matlabFunction this 1 x 6 row vector.

3. Second equation, differentiate using diff with respect to each state variable. That means use diff six times. Now, matlabFunction this 1 x 6 row vector.

4. Third equation, differentiate using diff with respect to each state variable. That means use diff six times. Now, matlabFunction this 1 x 6 row vector.

5. Fourth equation, differentiate using diff with respect to each state variable. That means use diff six times. Now, matlabFunction this 1 x 6 row vector, 1 entry at a time. That means use matlabFunction six times here.

6. Fifth equation, differentiate using diff with respect to each state variable. That means use diff six times. Now, matlabFunction this 1 x 6 row vector, 1 entry at a time. That means use matlabFunction six times here.

7. Sixth equation, differentiate using diff with respect to each state variable. That means use diff six times. Now, matlabFunction this 1 x 6 row vector, 1 entry at a time. That means use matlabFunction six times here.

That makes 21 numerical files.

In a new script file, put these all together to make the 6x6 Jacobian numerical function, and evaluate it numerically. Then, find eigenvalues.

I like this idea and think it might work.

Do you like this?

I think it’s cool.

[Noob]

Hi Walter!

I've currently split up my symbolic Jacobian function and stored them into 6 variables J1, J2, ... , J6, where each variable is a 1x6 row rector.

I'm also using expand( ), followed by simplify( ), on each of J1, ... , J6.

Then I'll matlabFunction each of J1, ... , J6.

I wonder if this could work.

Thanks!

[Walter Roberson]

"That makes 21 numerical files."

You have six equations, each differentiated by six variables, so you should end up with either 6 or 36 files, not 21.

[Noob]

Hi Walter!

I just used 21, 1 file each for the first 3 rows of the Jacobian -- these files are not complex at all.

Then 1 file for each entry of the last 3 rows of the Jacobian -- these files are probably still too complex.

When I set Optimize to false, all 21 files got written by matlabFunction in just a minute or two. Super, super fast.

However, I still got the message in those numerical files that, "the code analysis is too complex".

So now, I am using the default code optimization by matlabFunction.

If that doesn't work either, I might then try passing in "Sparse", true.

If everything fails, I'll learn to write a central difference method, maybe?

Or use that very popular one on the File Exchange, called Jacobest within the derivest file, I think.

If you have more thoughts, please feel free to share.

Thanks!

[Noob]

Hi Walter!

Here's the end of my code, if you wanted to take a look:

J = expand(J);
J = simplify(J);
J1 = J(1,:);
J2 = J(2,:);
J3 = J(3,:);
J41 = J(4,1);    J42 = J(4,2);    J43 = J(4,3);    J44 = J(4,4);  J45 = J(4,5);    J46 = J(4,6);   
J51 = J(5,1);    J52 = J(5,2);    J53 = J(5,3);    J54 = J(5,4);  J55 = J(5,5);    J56 = J(5,6); 
J61 = J(6,1);    J62 = J(6,2);    J63 = J(6,3);    J64 = J(6,4);  J65 = J(6,5);    J66 = J(6,6); 
%% create numerical Matlab functions and write them to separate files
matlabFunction(J1,'file','J1');
matlabFunction(J2,'file','J2');
matlabFunction(J3,'file','J3');
matlabFunction(J41,'file','J41',"Sparse",true);
matlabFunction(J42,'file','J42',"Sparse",true);
matlabFunction(J43,'file','J43',"Sparse",true);
matlabFunction(J44,'file','J44',"Sparse",true);
matlabFunction(J45,'file','J45',"Sparse",true);
matlabFunction(J46,'file','J46',"Sparse",true);
matlabFunction(J51,'file','J51',"Sparse",true);
matlabFunction(J52,'file','J52',"Sparse",true);
matlabFunction(J53,'file','J53',"Sparse",true);
matlabFunction(J54,'file','J54',"Sparse",true);
matlabFunction(J55,'file','J55',"Sparse",true);
matlabFunction(J56,'file','J56',"Sparse",true);
matlabFunction(J61,'file','J61',"Sparse",true);
matlabFunction(J62,'file','J62',"Sparse",true);
matlabFunction(J63,'file','J63',"Sparse",true);
matlabFunction(J64,'file','J64',"Sparse",true);
matlabFunction(J65,'file','J65',"Sparse",true);
matlabFunction(J66,'file','J66',"Sparse",true);

[Walter Roberson]

You should always specify the 'vars' parameter, as some of the expressions might happen to not involve all of the variables.

Your J41 and so on are scalar expressions; using 'Sparse', true is not going to help at all.

[Noob]

Hi Walter!

Could I maybe bypass matlabFunction altogether and work both symbolically and numerically?

That is, within my symbolic code, can I just input all parameter values, and state-variable values (say, using a fixed point I've found), using double-precision numbers such as a = 1.235592791?

Then use the symbolic Jacobian command. Presumably, this Jacobian is now considered "numerical", since all the inputs to it were double-precision numbers.

Then compute the eigenvalues of the numerical Jacobian, using eig( ).

Are there any issues with this workflow?

Thanks!

[Noob]

Hi Walter!

Seems like I should, at this point in time, use a finite difference method to compute a numerical Jacobian, and stop using symbolic math. I'll check out the highly popular Jacobest within the derivest file on the File Exchange soon.

Thanks again!

[Noob]

Hi Walter!

I got it to work!!!!!!!!!

By reducing the number of symbolic variables and parameters.

Just make all the parameters known to the symbolic code, specify them all, instead of later evaluating the Jacobian at so many variable values and parameter values.

matlabFunction worked fast, and the numerical file doesn't come with any warning message that the code is too complex to analyze.

If you see any issues with my specifying all the parameter values, please let me know.

Reducing the number of symoblic variables was the key here, I think.

[Noob]

Hi Walter!

Different from my Simulation workflow, for my Stability workflow, I should immediately specify all parameter values (likely, using double-precision numbers) in order to reduce the number of symbolic variables. This seems to have resolved my issues.

Thanks!

[Noob]

Hi Walter!

Just wanted to say thanks again! Today, I think I am deciding to compute the Jacobian function numerically, using a code from the File Exchange and / or write a central difference code myself. I'm afraid of making drastic changes to my workflow (which I've worked hard to debug over a long period of time) in order to accommodate a massive symbolic Jacobian.

Thanks again for your help!