Solve 2 equations with two unknowns and an array
2 ビュー (過去 30 日間)
古いコメントを表示
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)
0 件のコメント
回答 (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 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!