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

Question

Gabriel McQueen 2022-12-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1872617-solving-system-of-non-linear-equations-for-array-of-variables-gives-0x1-syms

评论： Walter Roberson 2022-12-6

在 MATLAB Online 中打开

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2022-12-6

在 MATLAB Online 中打开

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 2022-12-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1872617-solving-system-of-non-linear-equations-for-array-of-variables-gives-0x1-syms#answer_1121737

在 MATLAB Online 中打开

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