How to solve for a in a double nested loop.

3 次查看（过去 30 天）
sophp
sophp 2018-2-6
编辑： JK 2018-2-8
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.

请先登录，再回答此问题。

回答（2 个）

JK
JK 2018-2-7
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.
  1 个评论
sophp
sophp 2018-2-7
编辑：sophp 2018-2-7
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.

请先登录，再进行评论。

JK
JK 2018-2-8
编辑：JK 2018-2-8
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

类别

在 Help Center 和 File Exchange 中查找有关 Manual Performance Optimization 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by