How to solve 3 non-linear equations with 3 unknowns?

Good day everyone

I want to solve 3 non-linear equations with 3 unknown variables. I have written a function for this. When execute the function in the command window I run into this warning:

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.

> In sym/solve (line 304)

In myFunction (line 26)

ans =

[x*(exp(24424673397975836262400/(252133388810140916979*y)) - 1) - 1748446772768461537/225179981368524800 == 0, z + (28014820978904546331*y*exp(-24424673397975836262400/(252133388810140916979*y)))/(15144266739816366080*x) - 5949100062310601/9007199254740992 == 0, (721279627821056*z)/17034449605697695 + x*(exp((252133388810140916979*((34107876141310053*z)/25961484292674138142652481646100480 + 152237593508774139/6490371073168534535663120411525120))/(850705917302346158658436518579420528640*y)) - 1) + 6011786680494931611/340688992113953900 == 0]

Could anyone please assist me in what I am doing wrong? I need 1 numerical awnser for x,y,z.

Thanks

function F = myFunction (z)


q = 1.60218e-19         %electron discharge
V_oc = 4*44.8           %Open circuit voltage for 4 panels
N_s = 72                %number of cells connected in series
k = 1.38065e-23         %Boltzman constant
T = 298.15              %cell temperature
I_sc = 8.69             %short circuit current
R_sh = 193.6592708      %shunt resistance
dVdI_oc = 1/-1.514044   %Inverse of slope ate open circuit
V_mp = 4*36.6           %Maximum power voltage
I_mp = 8.2              %Maximum power current
I_pv = I_sc

%I_o = x
%n = y
%R_s = z

syms x y z

F(1) = ((x*(exp((q*V_oc)/(y*N_s*k*T))-1))-I_sc+(V_oc/R_sh))==0;
F(2) = (z + (dVdI_oc) + (1/((q*x/(y*N_s*k*T))*exp(q*V_oc/(y*N_s*k*T)))))==0;
F(3) = (x*((exp((q*(V_mp+(I_mp*z)))/y*N_s*k*T))-1)+(V_mp+(I_mp*z))/R_sh + I_mp + I_sc)==0;

sol = solve([F(1),F(2),F(3)],[x,y,z]);
end

[Alex Sha]

Hi, I get the results as below:

x: -7.86254507709137
y: -22.0861082560496
z: -416.744522416098

Answer 1

[Sam Chak]

Hi @Stefan Krieg,

You can use vpasolve().

q       = 1.60218e-19;      %electron discharge
V_oc    = 4*44.8;           %Open circuit voltage for 4 panels
N_s     = 72;               %number of cells connected in series
k       = 1.38065e-23;      %Boltzman constant
T       = 298.15;           %cell temperature
I_sc    = 8.69;             %short circuit current
R_sh    = 193.6592708;      %shunt resistance
dVdI_oc = 1/-1.514044;      %Inverse of slope ate open circuit
V_mp    = 4*36.6;           %Maximum power voltage
I_mp    = 8.2;              %Maximum power current
I_pv    = I_sc;
% I_o = x
% n = y
% R_s = z
syms x y z
F(1) = ((x*(exp((q*V_oc)/(y*N_s*k*T))-1))-I_sc+(V_oc/R_sh))==0;
F(2) = (z + (dVdI_oc) + (1/((q*x/(y*N_s*k*T))*exp(q*V_oc/(y*N_s*k*T)))))==0;
F(3) = (x*((exp((q*(V_mp+(I_mp*z)))/y*N_s*k*T))-1)+(V_mp+(I_mp*z))/R_sh + I_mp + I_sc)==0;
sol = vpasolve([F(1),F(2),F(3)],[x,y,z])

sol = struct with fields:
    x: -2.636463254746580127870476337846e+49
    y: 7.2714991697885760498634526499995e+48
    z: -416.7445224177292689999788029582

[Stefan Krieg]

Thank you, much appreciated. This is not the expected solution, I will have look at my equations again.

Answer 2

[Walter Roberson]

The numeric solutions being produced are dubious.

Notice how if we increase the number of digits to process at, that the magnitudes of x and y go up a lot, but z stays much the same. Then when we back-substitute, the remainder for the first equation stays much the same, but the remainder for the other two equations changes notably.

Can we do better?

Q = @(v) sym(v);   %convert to symbolic number

q       = Q(1.60218e-19);      %electron discharge
V_oc    = 4*Q(44.8);           %Open circuit voltage for 4 panels
N_s     = Q(72);               %number of cells connected in series
k       = Q(1.38065e-23);      %Boltzman constant
T       = Q(298.15);           %cell temperature
I_sc    = Q(8.69);             %short circuit current
R_sh    = Q(193.6592708);      %shunt resistance
dVdI_oc = 1/Q(-1.514044);      %Inverse of slope ate open circuit
V_mp    = 4*Q(36.6);           %Maximum power voltage
I_mp    = Q(8.2);              %Maximum power current
I_pv    = I_sc;
% I_o = x
% n = y
% R_s = z
syms x y z
F(1,1) = ((x*(exp((q*V_oc)/(y*N_s*k*T))-1))-I_sc+(V_oc/R_sh))==0;
F(2,1) = (z + (dVdI_oc) + (1/((q*x/(y*N_s*k*T))*exp(q*V_oc/(y*N_s*k*T)))))==0;
F(3,1) = (x*((exp((q*(V_mp+(I_mp*z)))/y*N_s*k*T))-1)+(V_mp+(I_mp*z))/R_sh + I_mp + I_sc)==0;

F

F = 

%cross-check 50
digits(50)
sol = vpasolve(F,[x,y,z])

sol = struct with fields:
    x: 5.778320712828577486739426085443008797321206724901e+66
    y: -1.5936901145812319849095567081810637894296223142398e+66
    z: -416.74452278539695297294985201961716638096566722981

subs(F, sol)

ans = 

%cross-check 100
digits(100)
sol = vpasolve(F,[x,y,z])

sol = struct with fields:
    x: 1.778606202685051504194240299371941120912764206302539905506042308878578946691993020069557647352316783e+117
    y: -4.905485732079848187309318642813357222474541837281366191242543485143184989714742976625190643174205242e+116
    z: -416.747365580410707307368431863945702771674846388716182429050639941886945690525098549691804405869203

subs(F, sol)

ans = 

%cross-check 200
digits(200)

sol = vpasolve(F,[x,y,z])

sol = struct with fields:
    x: 2.1253277269192755215287267224132413071758900048346594765947120534469673714767739880357567466412505581730887441263769539733683919983979248358548464957150821112414062811924734376764374869229717715842355e+217
    y: -5.7914206399459399153919106726969579062403404550768843089398315107147382482227843273334578765531375238518851528617878820724697219915991584359971210064143088616577428137964903429226049413464891654291869e+216
    z: -1068.1081891420746454645013840766059447546928609541615260330634049789077827490389047336826713843612196750876654685256342672965304555594108001981875797439687325925188950833936463137625828173999255179196

%cross-check
subs(F, sol)

ans = 

%can we do better? What happens if we solve iteratively?

y_partial = solve(F(1), y)

y_partial = 

F2 = subs(F(2:end), y, y_partial)

F2 = 

z_partial = solve(F2(1), z)

z_partial = 

F3 = simplify(subs(F2(2:end), z, z_partial))

F3 = 

xsol = vpasolve(F3, x)

xsol = 

zsol = subs(z_partial, x, xsol)

zsol = 

ysol = subs(y_partial, x, xsol)

ysol = 

subs(F, [x, y, z], [xsol, ysol, zsol])

ans = 

%that's pretty good!