How to solve for a in a double nested loop.

This is just a general question as I am new to MATLAB and would like to see how more experienced users would tackle this.

Say if I had 3 variables a, b and c. I wish to find the point where 0=g(a,b,c) for a between limits a_1 and a_2 and b within limits b_1 and b_2. For each value of a_x and b_x, I wish to solve 0=g(a,b,c) for the values of c_x (will be two solutions). How would you structure the solution to this problem? What functions would you use? I have tried via "solver" but it is too slow and does not find more than one solution.

Answer 1

[JK]

I think you must provide more detail to get a good answer. Is g a vector of three dimensions (a,b,c) or is the result a scalar? is it g(x,a,b,c)? are you looking for the parameters a, b, c or for x?

If you want to use bounds for a, b, and c, try least square solvers such as lsqnonlin. It will require the optimization toolbox.

To speed up the solver, calculate the jacobian matrices. It can be specified as a second input argument in the function handle. Set the options of fsolve to:

options = optimoptions('fsolve','SpecifyObjectiveGradient','on')

Matlab can calculate the Jacobian, if you have the symbolic toolbox installed.

to calculate x for known a, b, and c use fzero, which is faster than fsolve but only works for scalars in the input and output.

Matlab will only find one solution, you have to chose the initial point carefully. Try to simplify your equation by eliminating non-linear terms with minor influence to calculate a good starting point for your solver routine.

It will help if you could post your equation system.

[sophp]

To be more specific this is the code I have produced so far

%defining constants
V=8.5;
P=3e5; %Pa
R=8.3145; %kJ/mol.K
Cp=1.09 %kJ/kgK
Mr_air=29e-3; %kg/mol
dH_1=-1850; %kJ/mol
dH_2=-1423;
dH_3=-3273;
F_ao=0.1; %molar flow rate of benzene
Ti=600:1:850; 
y_A=0.001:0.0001:0.01;
hold on;
for i=1:numel(Ti)
      for j=1:numel(y_A)
          To=To       
          k1(i) = 1e7.*exp(-12700./To);
          k2(i) = 5e4.*exp(-10800./To);
          k3(i) = 7e7.*exp(-15000./To); 
          t(j)=(V.*y_A(j).*P)./(F_ao.*R.*To);
          X(i,j)=(t(j).*(k1(i)+k3(i)))./(1+t(j).*(k1(i)+k3(i)));
          Y_B(i,j) = (k1(i).*t(j))./(((k2(i).*t(j))+1).*(1+(t(j).*(k1(i)+k3(i)))))
          dH_gen(i,j)=-((V.*y_A(j)/F_ao.*P)./(R.*To)).*((((k1(i).*dH_1)+(k3(i).*dH_3)).*(1-X(j)))+(k2(i).*Y_B(j).*dH_2)); %enthalphy generated by exothermic reaction
          dH_rem(i,j)=(Mr_air./y_A(j)).*Cp.*(To-Ti(i));%enthalphy change of process fluid
          eqn(i,j) = dH_gen(i,j)-dH_rem(i,j);
          solx = fzero(eqn(i,j),To)
      end
  end

I am looping around Ti and y_A and using the functions defined above to find the value of To. Do I set To to a symbolic array? There may also be several values of To for a given value of Ti and y_A.

Answer 2

[JK]

sorry for the late answer. Here is the code:

%defining constants
tic

Ti=600:1:850; 
y_A=0.001:0.0001:0.01;

solx=zeros(numel(Ti),numel(y_A));

for i=1:numel(Ti)
      for j=1:numel(y_A)
          solx(i,j) = fzero(@(To) eqn(i,j,To),300);  % change 300 for 600 +i for differt solusions
      end
end

toc

figure
contour3(y_A,Ti,solx,50); 
xlabel('y_a');
ylabel('Ti');
zlabel('To');

function eq=eqn(i,j,To)  
    % Computation of Ti an yA from i an j
    Ti=600+i;
    y_A=0.001+0.0001*j;         % I am sure this can be done more elegant

      % Constants
      V=8.5;
      P=3e5; %Pa
      R=8.3145; %kJ/mol.K
      Cp=1.09; %kJ/kgK
      Mr_air=29e-3; %kg/mol
      dH_1=-1850; %kJ/mol
      dH_2=-1423;
      dH_3=-3273;
      F_ao=0.1; %molar flow rate of benzene

      k1 = 1e7.*exp(-12700./To);          
      k2 = 5e4.*exp(-10800./To);
      k3 = 7e7.*exp(-15000./To); 
      t=(V.*y_A.*P)./(F_ao.*R.*To);
      X=(t.*(k1+k3))./(1+t.*(k1+k3));
      Y_B = (k1.*t)./(((k2.*t)+1).*(1+(t.*(k1+k3))));
      dH_gen=-((V.*y_A/F_ao.*P)./(R.*To)).*((((k1.*dH_1)+(k3.*dH_3)).*(1-X))+(k2.*Y_B.*dH_2)); %enthalphy generated by exothermic reaction
      dH_rem=(Mr_air./y_A).*Cp.*(To-Ti);%enthalphy change of process fluid
      eq= dH_gen-dH_rem;
  end

it takes 8 sec to compute - even without jacobians.

I think this is an autothermic reactor you are modelling and it can have several operation temperautres. I used Ti as a starting point for the calculation of To. Change the starting point to get different temperatures. Make the system transient to examine how the temperautre will change over time and what happens if you chnage the inlet temperautre.

Compare the code to your original code, don't hesitate to ask further questions.

cheers

JK