solving system of non-linear equations for array of variables gives 0x1 syms

The code is supposed to output numbers for E1, E2, E3 but gives:

E =

struct with fields:

E1: [0×1 sym]

E2: [0×1 sym]

E3: [0×1 sym]

The third equations for E3 is a dependent equation based on the first two (eqn1-eqn2) if that may be a source of the problem.

Code is as follows:

Ep(0.05)

KeqPstandard = 

E = struct with fields:
    E1: [0×1 sym]
    E2: [0×1 sym]
    E3: [0×1 sym]

ans = struct with fields:
    E1: [0×1 sym]
    E2: [0×1 sym]
    E3: [0×1 sym]

function E = Ep(P)
    %Calculatation of extent of reaction
    epsilon = sym('E', [1 3]); %Setting extent of reaction as a variable to solve

    N0 = [1; %CO2 %Initial Moles of Each Species
          3; %H2
          0; %Me
          0; %H2O
          1];%CO

    N = [N0(1) - epsilon(1) - epsilon(2); %CO2 %Moles during reation of each species
         N0(2) - 3.*epsilon(1) - epsilon(2) - 2.*epsilon(3); %H2
         N0(3) + epsilon(1) + epsilon(3); %Me
         N0(4) + epsilon(1) + epsilon(3); %H2O
         N0(5) + epsilon(2) - epsilon(3)]; %CO

    NT = sum(N0) - 2*epsilon(1) - 2*epsilon(3); %Total moles

    y = [N/NT]; %molar ratio of each component

    KeqP=[(y(3) .* y(4)) ./ (y(1) * y(2) .^3 .*(P^2));
          (y(4) .* y(5)) ./ (y(1) .* y(2));
           y(3) ./ (y(5) .*y(2) .^2 .*P)];%Calculating Keq=product(yiP)^new(i)

    %Calculation of Keq based on P

    R = 8.314/1000; %kJ/(molK)
    T = 373;
    G = [3.0995.*log(P) - 23.038; %co2
         0; %h2
         3.0864.*log(P) - 3.8754; %ch3oh
         3.0971.*log(P) - 1.3816; %h2o
         3.1013.*log(P) - 29.689;]; %co


    deltaGrxn = [-G(1)-3.*G(2)+G(3)+G(4); %-CO2-3H2+Me+H2O
                -G(1)-G(2)+G(4)+G(5); %-Co2-H2+H2O+CO
                -G(5)-2.*G(2)+G(3)]; %-CO-2H2O+Me 
    
    KeqT=[exp(-deltaGrxn(1)/(R*T));
         exp(-deltaGrxn(2)/(R*T));
         exp(-deltaGrxn(3)/(R*T))];

    KeqPstandard = 0 == KeqP - KeqT%Relationship of thermodynamic and kinetic representations of Keq

    E = solve(KeqPstandard, epsilon)

end

[Walter Roberson]

KeqPstandard = 0 == KeqP - KeqT

if you stop the code right after that line, and you do

sol = vpasolve(vpa(rhs(KeqPstandard)))

then you get a series of 24 solutions. But if you substitute the solutions back into KeqPstandard, they turn out not to be good solutions to the real problem even if they might be acceptable to the approximated problem.

If you used vpasolve(rhs(KeqPstandard)) then you get a single solution that is [0 0.3 0] -- which gives a division by 0 when you substitute it back.

Answer 1

[Torsten]

The numerical solver does not find a solution.

Are you sure a solution exists ?

epsilon0 = rand(3,1);
P = 0.05;
options = optimset('MaxFunEvals',1000000,'MaxIter',1000000);
epsilon = fsolve(@(epsilon)Ep(epsilon,P),epsilon0,options);

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.000000e+06.

norm(Ep(epsilon,P))

ans = 184.6611

epsilon=epsilon.^2

epsilon = 3×1

0.0000 1.0000 0.0008

function KeqPstandard = Ep(epsilon,P)
epsilon = epsilon.^2;
    %Calculatation of extent of reaction
    %epsilon = sym('E', [1 3]); %Setting extent of reaction as a variable to solve
    N0 = [1; %CO2 %Initial Moles of Each Species
          3; %H2
          0; %Me
          0; %H2O
          1];%CO
    N = [N0(1) - epsilon(1) - epsilon(2); %CO2 %Moles during reation of each species
         N0(2) - 3.*epsilon(1) - epsilon(2) - 2.*epsilon(3); %H2
         N0(3) + epsilon(1) + epsilon(3); %Me
         N0(4) + epsilon(1) + epsilon(3); %H2O
         N0(5) + epsilon(2) - epsilon(3)]; %CO
    NT = sum(N0) - 2*epsilon(1) - 2*epsilon(3); %Total moles
    y = [N/NT]; %molar ratio of each component
    KeqP=[(y(3) .* y(4)) ./ (y(1) * y(2) .^3 .*(P^2));
          (y(4) .* y(5)) ./ (y(1) .* y(2));
           y(3) ./ (y(5) .*y(2) .^2 .*P)];%Calculating Keq=product(yiP)^new(i)
    %Calculation of Keq based on P
    R = 8.314/1000; %kJ/(molK)
    T = 373;
    G = [3.0995.*log(P) - 23.038; %co2
         0; %h2
         3.0864.*log(P) - 3.8754; %ch3oh
         3.0971.*log(P) - 1.3816; %h2o
         3.1013.*log(P) - 29.689;]; %co
    deltaGrxn = [-G(1)-3.*G(2)+G(3)+G(4); %-CO2-3H2+Me+H2O
                -G(1)-G(2)+G(4)+G(5); %-Co2-H2+H2O+CO
                -G(5)-2.*G(2)+G(3)]; %-CO-2H2O+Me 
    
    KeqT=[exp(-deltaGrxn(1)/(R*T));
         exp(-deltaGrxn(2)/(R*T));
         exp(-deltaGrxn(3)/(R*T))];
    KeqPstandard = KeqP - KeqT;%Relationship of thermodynamic and kinetic representations of Keq
end