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)

Answer 1

[hosein Javan]

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

[Vito Di Bernardo]

Hi Hosein Javan,

First of all, thank you for your answer!

But how should I implement the two equations and solve the two unknowns?

Thanks in advance for your answer!

[hosein Javan]

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]

Hi Hosein Javan,

I have now solved the problem by making a while loop, since I also had to implement a condition on the final outcome. But when I increment j by 1000 (to simulate faster), j as well as my final value V result in a matrix with 1000 colums all equal to zero except the first, which contains the correct value. The plot is now a line from every correct point to the origin, instead of the points being connected. My computer is now running the +1 code, I think this will be correct, but it takes a while on my fairly old machine.

Have you any idea what the problem might be?

Kind regards,

Vito

[hosein Javan]

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]

Hi,

The answer was correct and the plot as well (after several hours). I am a little concerned for posting my code online. I am working on my thesis and the information is quite confidential.

But I wanted to thank you again for the help you gave!

Cheers,

Vito

[hosein Javan]

don't mention it. best wishes.

[Vito Di Bernardo]

Thank you su much!