expression is too large, symengine can't calculate

I use symbolic expression to calculate sensitivity equasions of ODEs with 19 ODEs and 42 variables. And the error message is "the expression is too large , symengin can't calculate " when using it to calculate the second order sensitivity equations. I checked the dimension is more than 3 ten thousand and it isn't the sparse matrix. I debug into the symbolic there is no sourcecode, how do I do next?

Part codes like next:

     function obj = compute_2nd_order_sensitivity_equations(obj)
            M = obj.M;
            N = obj.N;
            if isempty(obj.S) || isempty(obj.dS) || isempty(obj.S0)
                obj = obj.compute_1st_order_sensitivity_equations;
            end  
            [SS,dSS,SS0] = sensitivity_equations_another([obj.x;obj.S(:)], obj.theta, [obj.dx;obj.dS(:)], [obj.x0;obj.S0(:)]); 
            SS(1:M,:)=[];  dSS(1:M,:)=[]; SS0(1:M,:)=[]; % remove the first order sensitivity equations
            obj.SS = reshape(SS,M,N,N) ;                   % reshape the 2nd order equations 
            obj.dSS = reshape(dSS,M,N,N);
            obj.SS0 = reshape(SS0,M,N,N);
            obj.ODE_2nd_order_sensitivity_rhs  = matlabFunction(obj.dx,obj.dS(:),obj.dSS(:),'vars', {'t', obj.x,obj.S(:),obj.SS(:), obj.theta});    % define the right hand side of the ODE as a function of time, system states + sensitivity, and parameter
            obj.ODE_2nd_order_sensitivity_init = matlabFunction(obj.x0,obj.S0(:),obj.SS0(:),'vars', {obj.theta});                                     % define the initial condition as a function of parameters
           end

symengine give error messages in function of matlabFunction

......

else

body = mup2matcell(funs, opts.Sparse);

body = renameInputs(body,vars,inputs);

g = symengine('makeFhandle',varnames,body);

end

[Walter Roberson]

It would be more common to use odeFunction() instead of matlabFunction in constructing equations ?

Answer 1

[Walter Roberson]

Create

eqns = [obj.dx,obj.dS(:);obj.dSS(:)];

Now break that up into chunks. For example:

At_a_time = 10;
num_eqn = length(eqns);
num_chunks = ceil(num_eqn ./ At_a_time);
fh = cell(num_chunks, 1);
for K = 1 : At_a_time : num_eqn - At_a_time
    fn{K} = matlabFunction(eqns(K:K+At_a_time-1), 'vars', {'t', obj.x,obj.S(:),obj.SS(:), obj.theta});
end
if num_chunks * At_a_time ~= num_eqn    %you should double-check the boundary conditions here
    fn{K+1} = matlabFunction(eqns(K+At_a_time:end), 'vars', {'t', obj.x,obj.S(:),obj.SS(:), obj.theta});
end

And now:

obj.ODE_2nd_order_sensitivity_rhs = @(varargin) cell2mat(cellfun(@(FH) FH(varargin{:}), fn, 'uniform', 0));

which should execute the handles one at a time and form a vector of the results.

Yes, this is a bit of a hack, but if it gets the job done then it gets the job done.

[Wang Hong]

@ Walter Roberson: Thanks a lot! I will try it to check whether it can continue to calculate

[Wang Hong]

@Walter Roberson, Though the origianl problem can be executed, the next steps face another obtacle. Owing to the result of is so different, when we compute ode45 with the function handle,it gives the error information is : the number of input parameters is not enough , and we call it like this:

[t,y] = ode45(@(t,x_values)obj.ODE_2nd_order_sensitivity_rhs(t,x_values,para),obj.time_points_to_evaluate,obj.ODE_2nd_order_sensitivity_init(para));

[Walter Roberson]

You use matlabFunction with 'vars' and provide a cell array with five elements. That is going to convert into a function that expects five parameters. However you are only passing three parameters to it. It needs t, x S, SS, and theta parameters.

Answer 2

[Wang Hong]

@Walter Roberson,

I tried different methods, such as:

1. p_values={obj.x,obj.S(:),obj.SS(:)};

options = odeset('RelTol',1e-1);

[t,y] = ode45(@(t,x_values)obj.ODE_2nd_order_sensitivity_rhs(t,obj.dx,obj.dS(:),obj.dSS(:),obj.para),obj.time_points_to_evaluate,obj.ODE_2nd_order_sensitivity_init(para),options,p_values);

and

2. p_values={obj.x,obj.S(:),obj.SS(:)};

options = odeset('RelTol',1e-1);

[t,y] = ode45(@(t,x_values)obj.ODE_2nd_order_sensitivity_rhs(t,x_values,obj.para),obj.time_points_to_evaluate,obj.ODE_2nd_order_sensitivity_init(para),options,p_values);

but the error information are both : too many input arguments, what makes me confusing are the bold parts ODE_2nd_order_sensitivity_rhs(t,obj.dx,obj.dS(:),obj.dSS(:),obj.para) and p_values, I think the p_values is the input arguments and what is the previous?

My God , I don't know how to pass the parameters

[Walter Roberson]

MATLAB is going to detect that you did not pass a nonlinear constraint function after the boundary conditions. It is then going to see the p_values as an extra parameter that it must pass to the function along with the regular two parameters, so it would pass it along with t and x_values to your anonymous function that you defined as only expecting two parameters.

The function you call does not appear to expect p_values as a parameter.

[Wang Hong]

Oh my God, I find the right function call form, but it shows another error: out of memory.

I tried two methods when I created the sensitivity equations, and I find the first require more memory

but it is only 53482688bit and far less than my free memory capacity.

1.obj.ODE_2nd_order_sensitivity_rhs = matlabFunction(obj.dx,obj.dS(:),obj.dSS(:),'vars', {'t', obj.x,obj.S(:),obj.SS(:), obj.theta});

2.obj.ODE_2nd_order_sensitivity_rhs = matlabFunction(obj.dx,obj.dS(:),obj.dSS(:),'vars', {'t', [obj.x;obj.S(:);obj.SS(:)], obj.theta});

Whether the ODEs sensitivity shouldn't be solved by symbolic?

[Walter Roberson]

The second form is more likely correct. However if you are constructing a function to pass to the ode* numeric routines then you should follow the flow shown in the first example for odeFunction()

[Wang Hong]

The two form corresponding different define form can be executed for small-scale ODEs . But for larg scale ODEs will throw different error information which makes it not to continue. The first tells me:out of memory and the second tell me :the symbolic expression is too large to calculate .

here the size of obj.x is 19*1; the size of obj.S() is 798*1 and the sizie of obj.SS(:) is 33516*1

and the original ODEs includes 19 equations and 42 variables. In ordre to calculate the geodesic path integration, we need fist order and second order sensitivity. Owing to the implicit style of ODEs, we need to calculate the symbolic expression of the right side

[Walter Roberson]

When you list SS as part of var option you are declaring that the system has over 33 thousand variables. You cannot conduct a sensitivity analysis of a system with over 33 thousand variables in any reasonable amount of time and space. Not even if those 33 thousand are to be considered constants for the purpose of the equations. For a sensitivity analysis you would be wanting to create a Jacobian so that you can analyze how changes in two variables at a time influence the equations, and that gets you to a 34000-ish by 34000-ish symbolic array that would be potentially quite complicated.

[Wang Hong]

Thank you for your suggestion. I will try it! But I am curious about how to do sensitivity analysis when the number of variable is large

[Walter Roberson]

You don't do a sensitivity analysis when the number of variables is quite large. Not unless you are only concerned about the change with respect to one variable at a time, not about sensitivity as multiple variables change.

If you only care about one variable at a time then differentiate each equation with respect to each variable. I know that there is a function to do that but the function name is not coming to mind. This would not involve the use of matlabFunction. matlabFunction would only be for numeric work but you want analytic work.

Possibly you would be willing to establish a base numeric point and analyze the effect of changing one variable at a time at that particular numeric point. That would involve looping substituting all the numeric values except that one variable at a time you substitute numeric value plus symbolic delta. In some cases you will get a polynomial output; others might be difficult to analyze, such as if the equations involved the bessel functions.