Solve 2 equations with two unknowns and an array
显示 更早的评论
Should I use a for-loop instead of an array(preferred)? (2018a)
First try (When I put j equal to a roundom value between 1 and 1000, matlab yields a solution):
j = 1:1:1000;
T = 343;
n_H2O = 2;
alpha = 0.5;
j_0 = 1;
P_a = 3;
P_c = 3;
P_SAT = 0.307;
X_o2_d = 0.19;
X_H2O_a = 0.1;
X_H2O_d = 0.1;
t_a = 0.00035;
t_c = 0.00035;
t_M = 0.000125;
D_eff_O2_H2O = 0.000149;
D_eff_H2_H2O = 0.0000295;
F = 96485;
R = 8.314;
D_lambda = 0.00000381;
n_SAT_drag = 2.5;
M_m = 1;
rho_dry = 0.00197;
syms x
syms y
a_w_anode = (P_a./P_SAT).*(X_H2O_a-((t_a.*alpha_unknown.*j.*R.*T)./(n_H2O.*F.*P_a.*101325.*D_eff_H2_H2O)));
a_w_cathode = (P_c./P_SAT).*(X_H2O_d-((t_c.*(1+alpha_unknown).*j.*R.*T)./(n_H2O.*F.*P_c.*101325.*D_eff_O2_H2O)));
lambda_anode_1 = 14.*a_w_anode;
lambda_cathode_1 = 10+4.*a_w_cathode;
lambda_anode_2 = ((11.*alpha_unknown)./n_SAT_drag)+C;
lambda_cathode_2 = ((11.*alpha_unknown)./n_SAT_drag)+C.*exp((j.*M_m.*n_SAT_drag.*t_M)./(22.*F.*rho_dry.*D_lambda));
Eq1 = 14.*a_w_anode == ((11.*alpha_unknown)./n_SAT_drag)+C;
Eq2 = 10+4.*a_w_cathode == ((11.*alpha_unknown)./n_SAT_drag)+C.*exp((j.*M_m.*n_SAT_drag.*t_M)./(22.*F.*P_c.*D_lambda))
[x,y] = solve(Eq1,Eq2,[alpha_unknown,C])
alpha_unknown_sol = double(x)
C_sol = double(y)
%lambda(z) = ((11.*alpha_unknown_sol)./n_SAT_drag)+C_sol.*exp((j.*M_m.*n_SAT_drag.*z)./(22.*F.*P_c.*D_lambda));
%sigma(Z) = (0.005193.*lambda(z)-0.00326).*exp(1268.*(1./303-1/T));
fun = @(z) 1./(0.005193.*((11.*alpha_unknown_sol)./n_SAT_drag)+C_sol.*exp((j.*M_m.*n_SAT_drag.*z)./(22.*F.*P_c.*D_lambda))-0.00326).*exp(1268.*(1./303-1./T))
ASR_m = integral(fun,0,t_M,'ArrayValued',true)
V_ohmic = j.*ASR_m
plot(j,V_ohmic,'-*');
second try:
j = 1:1:1000;
T = 343;
n_H2O = 2;
alpha = 0.5;
j_0 = 1;
P_a = 3;
P_c = 3;
P_SAT = 0.307;
X_o2_d = 0.19;
X_H2O_a = 0.1;
X_H2O_d = 0.1;
t_a = 0.00035;
t_c = 0.00035;
t_M = 0.000125;
D_eff_O2_H2O = 0.000149;
D_eff_H2_H2O = 0.0000295;
F = 96485;
R = 8.314;
D_lambda = 0.00000381;
n_SAT_drag = 2.5;
M_m = 1;
rho_dry = 0.00197;
syms x
syms y
a_w_anode = (P_a./P_SAT).*(X_H2O_a-((t_a.*alpha_unknown.*j.*R.*T)./(n_H2O.*F.*P_a.*101325.*D_eff_H2_H2O)));
a_w_cathode = (P_c./P_SAT).*(X_H2O_d-((t_c.*(1+alpha_unknown).*j.*R.*T)./(n_H2O.*F.*P_c.*101325.*D_eff_O2_H2O)));
lambda_anode_1 = 14.*a_w_anode;
lambda_cathode_1 = 10+4.*a_w_cathode;
lambda_anode_2 = ((11.*alpha_unknown)./n_SAT_drag)+C;
lambda_cathode_2 = ((11.*alpha_unknown)./n_SAT_drag)+C.*exp((j.*M_m.*n_SAT_drag.*t_M)./(22.*F.*rho_dry.*D_lambda));
syms t
Eq1 = 14.*a_w_anode == ((11.*alpha_unknown)./n_SAT_drag)+C;
Eq2 = 10+4.*a_w_cathode == ((11.*alpha_unknown)./n_SAT_drag)+C.*exp((t.*M_m.*n_SAT_drag.*t_M)./(22.*F.*P_c.*D_lambda))
sol = solve([Eq1,Eq2],[x,y])
x = subs(sol.x, t, j)
y = subs(sol.y, t, j)
alpha_unknown_sol = double(x)
C_sol = double(y)
%lambda(z) = ((11.*alpha_unknown_sol)./n_SAT_drag)+C_sol.*exp((j.*M_m.*n_SAT_drag.*z)./(22.*F.*P_c.*D_lambda));
%sigma(Z) = (0.005193.*lambda(z)-0.00326).*exp(1268.*(1./303-1/T));
fun = @(z) 1./(0.005193.*((11.*alpha_unknown_sol)./n_SAT_drag)+C_sol.*exp((j.*M_m.*n_SAT_drag.*z)./(22.*F.*P_c.*D_lambda))-0.00326).*exp(1268.*(1./303-1./T))
ASR_m = integral(fun,0,t_M,'ArrayValued',true)
V_ohmic = j.*ASR_m
plot(j,V_ohmic,'-*');
The error with first try:
Error using plot
Vectors must be the same length.
Error in Ohmic_losses (line 56)
plot(j,V_ohmic,'-*');
The error with second try:
Error using symengine
Number of elements in NEW must match number in OLD.
Error in sym/subs>mupadsubs (line 160)
G = mupadmex('symobj::fullsubs',F.s,X2,Y2);
Error in sym/subs (line 145)
G = mupadsubs(F,X,Y);
Error in Ohmic_losses (line 43)
x = subs(sol.x, t, j)
回答（1 个）
hosein Javan
2020-8-10
you are using symbolic toolbox. I recommend to make a function file looking like this:
function y=eq(x)
% enter equations
% the equations must be in the form of f(x)=0
% for example if sin(x)=2*x then write y(1)=sin(x)-2*x
end
then use the following to solve the system numerically:
x0 = [0,0]; % initial guess of unknowns x1 and x2
x_sol = fsolve(@eq,x0) % solutions
7 个评论
Vito Di Bernardo
2020-8-11
hosein Javan
2020-8-11
编辑：Walter Roberson
2020-8-11
hello and you're welcom. I'll give you an example:
suppose you want to solve this system:
x(1) = sin(x(2))
x(2) = - exp(x(1))
your code will be:
eq = @(x) [
x(1) - sin(x(2))
x(2) + exp(x(1))
];
x0 = [0,0];
x = fsolve(eq,x0)
now in command window a message would appear that the solution is converged:
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
x =
-0.5469 -0.5787
notice that initial guess is sometimes important for functions that have bad behaviour. give it a try, if you didn't get your answer, we would discuss about it.
Vito Di Bernardo
2020-8-11
hosein Javan
2020-8-11
hello again. I must see and check the equations to say anything. I think you have write an iteration-based solver code yourself. if you know that the answer is correct then the job is done. but if not use "fsolve". if the solution is not found easily, make a 3d surface plot to see for which values, the answer is near zero to make a better initial guess. if the problem persists, show me your equations.
Vito Di Bernardo
2020-8-12
hosein Javan
2020-8-12
don't mention it. best wishes.
Vito Di Bernardo
2020-8-15
类别
在 帮助中心 和 File Exchange 中查找有关 Code Performance 的更多信息
2020-8-10
2020-8-15
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
