Strange error message in for loop code

Hello All,

I am getting an array assignment error that does not make any sense to me. My code, shown below, outputs a 1x141 array of values for the five primary species of gas that constitute air at equilibrium.

The code solves a system of five equations for the five unknowns. I am using vpasolve because several of the equations are nonlinear.

I'm using a for loop to solve the eqns, which are temperature dependent, over a range of temperatures, and saving the values of the 5 unknowns in the 1x141 arrays so that they can be plotted versus temperature.

For some reason, my code runs perfectly fine up until the 121st iteration in my loop, then it gives me the following error:

In an assignment A(I) = B, the number of elements in B and I must be the same.

Error in Untitled (line 42) O2(n) = double((solnO2/Total_Moles)*100);

I don't understand why I am getting this error.

Can anyone shed any light on this?

Thanks very much.

clear all
close all
clc

syms nO nO2 nN nN2 nNO

C = [1.2e3, 18, 4];
eta = [-0.5, 0, 0];
theta = [59500, 113000, 75500];
rho = [10^(-3), 10^(-4), 10^(-5)];
nO2_0 = 0.00733;
nN2_0 = 0.0273;

Temp = (2000:100:16000);

O2 = zeros(1,141);
N2 = zeros(1,141);
O  = zeros(1,141);
N  = zeros(1,141);
NO = zeros(1,141);

n = 0;

for T = 2000:100:16000

      n = n + 1;

      K1 = C(1)*T^eta(1)*exp(-theta(1)/T);
      K2 = C(2)*T^eta(2)*exp(-theta(2)/T);
      K3 = C(3)*T^eta(3)*exp(-theta(3)/T);

      [solnO2, solnN2, solnO, solnN, solnNO] =... 
      vpasolve([nO == sqrt(nO2*(1/rho(1))*K1),...
      nN == sqrt(nN2*(1/rho(1))*K2),...
      nNO == (nN*nO)/((1/rho(1))*K3),...
       nO2 == 0.5*(2*nO2_0 - nO - nNO),...
       nN2 == 0.5*(2*nN2_0 - nN - nNO)],...
       [nO2, nN2, nO, nN, nNO]);

      Total_Moles = solnO2 + solnN2 + solnO + solnN + solnNO;

      O2(n) = double((solnO2/Total_Moles)*100);
      N2(n) = double((solnN2/Total_Moles)*100);
      O(n)  = double((solnO/Total_Moles)*100);
      N(n)  = double((solnN/Total_Moles)*100);
      NO(n) = double((solnNO/Total_Moles)*100);

end

Answer 1

[Roger Stafford]

I think what this means is that 'vpasolve' found two or more solutions on that 121st iteration, but it was trying to stuff these into a single position in O2(n). It would be best in your loop to check the size of your solution set and make the proper provision for more than one solution at a time:

   t = double((solnO2/Total_Moles)*100);
   nn = numel(tt);
   O2(n:n+nn-1) = t;
   ...
   n = n + nn; % Instead of n = n + 1

[Gabriel]

Thanks for your response.

In the second line of your code, did you mean to type

nn = numel(t);

Also, why would this vpasolve only start to find two or more solutions on the 121st iteration?

[Roger Stafford]

Yes I meant nn = numel(t);. Sorry.

As to why it only found multiple solutions at the 121st iteration, I can only say this. Your equations can be reduced to the solution to a sixth-degree polynomial. However, because of the presence of 'sqrt' in the original equations, not all six solutions to this polynomial will be valid solutions to those original equations. Presumably as you varied the parameter T, it caused K1, K2, and K3 to vary in such a way that at a certain point more of the polynomial's answers became valid than just one of them. I can't be more specific than that because I have not studied that polynomial in detail to see how its solutions relate to your original equations.

[Gabriel]

Ok, thanks for the suggestion.

If vpasolve gives me multiple solutions, and I am only interested in one, is there any reason why the following wouldn't work?

    Total_Moles = solnO2(1) + solnN2(1) + solnO(1) + solnN(1) + solnNO(1);

      O2(n) = double((solnO2(1)/Total_Moles));
      N2(n) = double((solnN2(1)/Total_Moles));
      O(n)  = double((solnO(1)/Total_Moles));
      N(n)  = double((solnN(1)/Total_Moles));
      NO(n) = double((solnNO(1)/Total_Moles));

I tried this, but am still getting the same array error.