Why do I get an error with the Solve function

I have a code where I am trying to performed constrained optimization using a Lagrangian. I have set up a series of first order conditions using the diff function. However when I attempt to solve the a system of 11 equations with 11 unknowns I get the following errors:

Error using mupadengine/feval (line 166) Out of memory.

Error in solve (line 293) sol = eng.feval('solve', eqns, vars, solveOptions);

I do not get an error if I comment out the solve function.

Any Advice? Code is Copied Below:

if true
  % code
timebudget=200;
initial1_200=.5;
initial1_400=.4;
initial1_600=.3;
initial1_800=.2;
initial1_1000=.1;
initial2_200=.55;
initial2_400=.45;
initial2_600=.35;
initial2_800=.25;
initial2_1000=.15;
initial3_200=.51;
initial3_400=.41;
initial3_600=.31;
initial3_800=.21;
initial3_1000=.11;
initial4_200=.45;
initial4_400=.4;
initial4_600=.35;
initial4_800=.3;
initial4_1000=.15;
initial5_200=.52;
initial5_400=.42;
initial5_600=.32;
initial5_800=.22;
initial5_1000=.12;
initial6_200=.6;
initial6_400=.5;
initial6_600=.4;
initial6_800=.3;
initial6_1000=.2;
initial7_200=.5;
initial7_400=.4;
initial7_600=.3;
initial7_800=.2;
initial7_1000=.1;
initial8_200=.53;
initial8_400=.43;
initial8_600=.33;
initial8_800=.23;
initial8_1000=.13;
initial9_200=.56;
initial9_400=.46;
initial9_600=.36;
initial9_800=.26;
initial9_1000=.16;
initial10_200=.8;
initial10_400=.58;
initial10_600=.38;
initial10_800=.3;
initial10_1000=.25;
r_200 = .1;
r_400 = .08;
r_600 = .06;
r_800 = .04;
r_1000 = .02;
syms t1;
syms t2;
syms t3;
syms t4;
syms t5;
syms t6;
syms t7;
syms t8;
syms t9;
syms t10;
syms lambda;
EE1 = (200/(1+((1-initial1_200)/initial1_200)*exp(-r_200*t1)))+(400/(1+((1-initial1_400)/initial1_400)*exp(-r_400*t1)))+(600/(1+((1-initial1_600)/initial1_600)*exp(-r_600*t1)))+(800/(1+((1-initial1_800)/initial1_800)*exp(-r_800*t1)))+(1000/(1+((1-initial1_1000)/initial1_1000)*exp(-r_1000*t1)));
EE2 = (200/(1+((1-initial2_200)/initial2_200)*exp(-r_200*t2)))+(400/(1+((1-initial2_400)/initial2_400)*exp(-r_400*t2)))+(600/(1+((1-initial2_600)/initial2_600)*exp(-r_600*t2)))+(800/(1+((1-initial2_800)/initial2_800)*exp(-r_800*t2)))+(1000/(1+((1-initial2_1000)/initial2_1000)*exp(-r_1000*t2)));
EE3 = (200/(1+((1-initial3_200)/initial3_200)*exp(-r_200*t3)))+(400/(1+((1-initial3_400)/initial3_400)*exp(-r_400*t3)))+(600/(1+((1-initial3_600)/initial3_600)*exp(-r_600*t3)))+(800/(1+((1-initial3_800)/initial3_800)*exp(-r_800*t3)))+(1000/(1+((1-initial3_1000)/initial3_1000)*exp(-r_1000*t3)));
EE4 = (200/(1+((1-initial4_200)/initial4_200)*exp(-r_200*t4)))+(400/(1+((1-initial4_400)/initial4_400)*exp(-r_400*t4)))+(600/(1+((1-initial4_600)/initial4_600)*exp(-r_600*t4)))+(800/(1+((1-initial4_800)/initial4_800)*exp(-r_800*t4)))+(1000/(1+((1-initial4_1000)/initial4_1000)*exp(-r_1000*t4)));
EE5 = (200/(1+((1-initial5_200)/initial5_200)*exp(-r_200*t5)))+(400/(1+((1-initial5_400)/initial5_400)*exp(-r_400*t5)))+(600/(1+((1-initial5_600)/initial5_600)*exp(-r_600*t5)))+(800/(1+((1-initial5_800)/initial5_800)*exp(-r_800*t5)))+(1000/(1+((1-initial5_1000)/initial5_1000)*exp(-r_1000*t5)));
EE6 = (200/(1+((1-initial6_200)/initial6_200)*exp(-r_200*t6)))+(400/(1+((1-initial6_400)/initial6_400)*exp(-r_400*t6)))+(600/(1+((1-initial6_600)/initial6_600)*exp(-r_600*t6)))+(800/(1+((1-initial6_800)/initial6_800)*exp(-r_800*t6)))+(1000/(1+((1-initial6_1000)/initial6_1000)*exp(-r_1000*t6)));
EE7 = (200/(1+((1-initial7_200)/initial7_200)*exp(-r_200*t7)))+(400/(1+((1-initial7_400)/initial7_400)*exp(-r_400*t7)))+(600/(1+((1-initial7_600)/initial7_600)*exp(-r_600*t7)))+(800/(1+((1-initial7_800)/initial7_800)*exp(-r_800*t7)))+(1000/(1+((1-initial7_1000)/initial7_1000)*exp(-r_1000*t7)));
EE8 = (200/(1+((1-initial8_200)/initial8_200)*exp(-r_200*t8)))+(400/(1+((1-initial8_400)/initial8_400)*exp(-r_400*t8)))+(600/(1+((1-initial8_600)/initial8_600)*exp(-r_600*t8)))+(800/(1+((1-initial8_800)/initial8_800)*exp(-r_800*t8)))+(1000/(1+((1-initial8_1000)/initial8_1000)*exp(-r_1000*t8)));
EE9 = (200/(1+((1-initial9_200)/initial9_200)*exp(-r_200*t9)))+(400/(1+((1-initial9_400)/initial9_400)*exp(-r_400*t9)))+(600/(1+((1-initial9_600)/initial9_600)*exp(-r_600*t9)))+(800/(1+((1-initial9_800)/initial9_800)*exp(-r_800*t9)))+(1000/(1+((1-initial9_1000)/initial9_1000)*exp(-r_1000*t9)));
EE10 = (200/(1+((1-initial10_200)/initial10_200)*exp(-r_200*t10)))+(400/(1+((1-initial10_400)/initial10_400)*exp(-r_400*t10)))+(600/(1+((1-initial10_600)/initial10_600)*exp(-r_600*t10)))+(800/(1+((1-initial10_800)/initial10_800)*exp(-r_800*t10)))+(1000/(1+((1-initial10_1000)/initial10_1000)*exp(-r_1000*t10)));
EE_Total= EE1 + EE2 +EE3+EE4+EE5+EE6+EE7+EE8+EE9+EE10
EE_Total_Constrained = EE_Total + lambda*(timebudget - t1 -t2 -t3 -t4 -t5 -t6 -t7 -t8 -t9 -t10)
FOC_1 = diff(EE_Total_Constrained, t1)==0;
FOC_2 = diff(EE_Total_Constrained, t2)==0;
FOC_3 = diff(EE_Total_Constrained, t3)==0;
FOC_4 = diff(EE_Total_Constrained, t4)==0;
FOC_5 = diff(EE_Total_Constrained, t5)==0;
FOC_6 = diff(EE_Total_Constrained, t6)==0;
FOC_7 = diff(EE_Total_Constrained, t7)==0;
FOC_8 = diff(EE_Total_Constrained, t8)==0;
FOC_9 = diff(EE_Total_Constrained, t9)==0;
FOC_10 = diff(EE_Total_Constrained, t10)==0;
FOC_lambda = diff(EE_Total_Constrained, lambda)==0;
solve([FOC_1, FOC_2, FOC_3, FOC_4, FOC_5, FOC_6, FOC_7, FOC_8, FOC_9, FOC_10, FOC_lambda], [t1, t2, t3, t4 ,t5, t6, t7, t8,t9, t10]);
% t1_Solve = sol.t1
% t2_Solve = sol.t2
% t3_Solve = sol.t3
% t4_Solve = sol.t4
% t5_Solve = sol.t5
% t6_Solve = sol.t6
% t7_Solve = sol.t7
% t8_Solve = sol.t8
% t9_Solve = sol.t9
% t10_Solve = sol.t10
% t1_Solve
% t2_Solve
% t3_Solve
% t4_Solve 
% t6_Solve 
% t7_Solve 
% t8_Solve
% t9_Solv
% t10_Solve
end

Answer 1

[Alan Weiss]

I suggest that you not try to solve this symbolically, but instead look for a numerical solution, using fsolve or some other way.

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 2

[John D'Errico]

Not every problem has a symbolic solution. In fact, most such problem will not. It is trivial to write problems that do not have one. And even if a symbolic solution exists on paper, it is also easy to create a problem where it is still theoretically impossible to derive that solution.

The problem is, MATLAB does not know that. In fact, you cannot know that no symbolic solution exists for some random, arbitrary problem. So MATLAB tries to find a solution, because there MIGHT be a solution to be found. It can spend a great deal of effort. And systems with MANY unknowns (11 is many in this context, for a symbolic problem) will create highly complex, messy terms in the search for that solution. Since problems with zillions of terms take a great deal of time to work with as well as a great deal of memory to store the intermediate results, it is very easy to blow out your available memory.

The solution, as Alan suggests, is to avoid a symbolic solver when that route is almost certain to fail. Use a tool like fsolve instead, a purely numerical solver.

Answer 3

[Manav]

function viewSolid(zvar, F, G, yvar, f, g, xvar, a, b)

%VIEWSOLID is a version for MATLAB of the routine on page 161

% of "Multivariable Calculus and Mathematica" for viewing the region

% bounded by two surfaces for the purpose of setting up triple integrals.

% The arguments are entered from the inside out.

% There are two forms of the command --- either f, g,

% F, and G can be vectorized functions, or else they can

% be symbolic expressions. xvar, yvar, and zvar can be

% either symbolic variables or strings.

% The variable xvar (x, for example) is on the

% OUTSIDE of the triple integral, and goes between CONSTANT limits a and b.

% The variable yvar goes in the MIDDLE of the triple integral, and goes

% between limits which must be expressions in one variable [xvar].

% The variable zvar goes in the INSIDE of the triple integral, and goes

% between limits which must be expressions in two

% variables [xvar and yvar]. The lower surface is plotted in red, the

% upper one in blue, and the "hatching" in cyan.

%

% Examples: viewSolid(z, 0, (x+y)/4, y, x/2, x, x, 1, 2)

% gives the picture on page 163 of "Multivariable Calculus and Mathematica"

% and the picture on page 164 of "Multivariable Calculus and Mathematica"

% can be produced by

% viewSolid(z, x^2+3*y^2, 4-y^2, y, -sqrt(4-x^2)/2, sqrt(4-x^2)/2, ...

% x, -2, 2,)

% One can also type viewSolid('z', @(x,y) 0, ...

% @(x,y)(x+y)/4, 'y', @(x) x/2, @(x) x, 'x', 1, 2)

%

if isa(f, 'sym') % case of symbolic input

ffun=inline(vectorize(f+0*xvar),char(xvar));

gfun=inline(vectorize(g+0*xvar),char(xvar));

Ffun=inline(vectorize(F+0*xvar),char(xvar),char(yvar));

Gfun=inline(vectorize(G+0*xvar),char(xvar),char(yvar));

oldviewSolid(char(xvar), double(a), double(b), ...

char(yvar), ffun, gfun, char(zvar), Ffun, Gfun)

else

oldviewSolid(char(xvar), double(a), double(b), ...

char(yvar), f, g, char(zvar), F, G)

end

%%%%%%% subfunction goes here %%%%%%

function oldviewSolid(xvar, a, b, yvar, f, g, zvar, F, G)

for counter=0:20

xx = a + (counter/20)*(b-a);

YY = f(xx)*ones(1, 21)+((g(xx)-f(xx))/20)*(0:20);

XX = xx*ones(1, 21);

%% The next lines inserted to make bounding curves thicker.

widthpar=0.5;

if counter==0, widthpar=2; end

if counter==20, widthpar=2; end

%% Plot curves of constant x on surface patches.

plot3(XX, YY, F(XX, YY).*ones(1,21), 'r', 'LineWidth', widthpar);

hold on

plot3(XX, YY, G(XX, YY).*ones(1,21), 'b', 'LineWidth', widthpar);

end;

%% Now do the same thing in the other direction.

XX = a*ones(1, 21)+((b-a)/20)*(0:20);

%% Normalize sizes of vectors.

YY=0:2; ZZ1=0:20; ZZ2=0:20;

for counter=0:20,

%% The next lines inserted to make bounding curves thicker.

widthpar=0.5;

if counter==0, widthpar=2; end

if counter==20, widthpar=2; end

for i=1:21,

YY(i)=f(XX(i))+(counter/20)*(g(XX(i))-f(XX(i)));

ZZ1(i)=F(XX(i),YY(i));

ZZ2(i)=G(XX(i),YY(i));

end;

plot3(XX, YY, ZZ1, 'r', 'LineWidth',widthpar);

plot3(XX, YY, ZZ2, 'b', 'LineWidth',widthpar);

end;

%% Now plot vertical lines.

for u = 0:0.2:1,

for v = 0:0.2:1,

x=a + (b-a)*u; y = f(a + (b-a)*u) +(g(a + (b-a)*u)-f(a + (b-a)*u))*v;

plot3([x, x], [y, y], [F(x,y), G(x, y)], 'c');

end;

end;

xlabel(xvar)

ylabel(yvar)

zlabel(zvar)

hold off