I am trying to solve the non-linear solver 'fsolve' to solve the rate equations. But the f(x) shows unprecedented high values and its shows no solution found.

Codes are as follows:

close all;
clear all;
clc;

% Define the function to be solved
fun = @(k) Rate_eqns(k);

% Set the initial guess
k0 = [1.6e18, 0, 0];

% Set the options for the solver
options = optimoptions('fsolve','Display','iter');



% Call the solver
[k,~,exitflag] = fsolve(fun,k0,options);

                                             Norm of      First-order   Trust-region
 Iteration  Func-count     ||f(x)||^2           step       optimality         radius
     0          4         5.16979e+58                        3.23e+40              1
     1          5         5.16979e+58              1         3.23e+40              1
     2          6         5.16979e+58           0.25         3.23e+40           0.25

No solution found.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the value of the function tolerance.

% Check the exit flag
if exitflag <= 0
    error('fsolve did not converge to a solution.')
end

fsolve did not converge to a solution.

% Display the results
fprintf('Nqw = %f\n', k(1));
fprintf('Nb = %f\n', k(2));
fprintf('Np = %f\n', k(3));

the function is:

function fcns = Rate_eqns(k)

Nqw = k(1);
Nb = k(2);
Np = k(3);

% Ninj = 2e18;
Isch = 15.8e-3; q = 1.6e-19; V = 1e-10;
Ninj = Isch/(q*V);
Ab = 1e8; Bb = 1.2e-10;
tau_b = 1/(Ab+Bb*Nb);
Aqw = 1e8; Bqw = 1.5e-10; Cqw = 5e-30;
tau_r = 1/(Bqw*Nqw); tau_nr = 1/(Aqw+Cqw*Nqw^2); 
tau_e = 10e-12;
tau_bw = 1e-12;
g = 285; alpha_i = 1.25; alpha_m = 16; ng = 3.42;
vg = 3e10/ng;
tau_p = 1/(vg*(alpha_i+alpha_m));


fcns(1) = Ninj - (Nb/tau_b)- (Nb/tau_bw) + (Nqw/tau_e);
fcns(2) = (Nb/tau_bw) -(Nqw/tau_r)- (Nqw/tau_nr)- (Nqw/tau_e) - vg*g*Np;
fcns(3) = vg*g - (Np/tau_p);



end

[Torsten]

Your code has no syntax errors.

So the problem is due to the equations you are trying to solve. So either your equations are wrong, the implementation of the equations is wrong or the solution of the equations is inherently difficult.

In all three cases, we are not able to help you because you are the only person who understands the physical background of what you try to do.

[Mohammad Fozlay]

Do you have any idea why ||f(x)||^2 is so large?

[Torsten]

Either because your parameters/equations/implementation are wrong or because your initial guess for Nqw, Nb and Np as [1.6e18, 0, 0] is far off.

You work with very small and very big numbers in your code - thus you can imagine that even if everything is set correctly, the problem will be very sensitive and difficult to solve.

Answer 1

[Torsten]

Here is a symbolic solution of your problem. It seems to boil down to determine the roots of a polynomial of order 6. Since you want positive solutions, you get

Nqw = 2152699271590929986.667971525434

Nb = 216230193480156287.19357257820151

Np = 16.521739130434782608695652173913

with the parameters and equations you wrote.

syms Nqw Nb Np Isch q V Ab Bb Aqw Bqw Cqw tau_e tau_bw g alpha_i alpha_m ng const
Ninj = Isch/(q*V);
tau_b = 1/(Ab+Bb*Nb);
tau_r = 1/(Bqw*Nqw); 
tau_nr = 1/(Aqw+Cqw*Nqw^2); 
vg = const/ng;
tau_p = 1/(vg*(alpha_i+alpha_m));
eqn1 = Ninj - (Nb/tau_b)- (Nb/tau_bw) + (Nqw/tau_e) == 0;
eqn2 = (Nb/tau_bw) -(Nqw/tau_r)- (Nqw/tau_nr)- (Nqw/tau_e) - vg*g*Np == 0;
eqn3 = vg*g - (Np/tau_p) == 0;
sol = solve([eqn1 eqn2 eqn3],[Nqw Nb Np])

sol = struct with fields:
    Nqw: [6×1 sym]
     Nb: [6×1 sym]
     Np: [6×1 sym]

Ischnum = 15.8e-3; 
qnum = 1.6e-19; 
Vnum = 1e-10;
Abnum = 1e8; 
Bbnum = 1.2e-10;
Aqwnum = 1e8; 
Bqwnum = 1.5e-10; 
Cqwnum = 5e-30;
tau_enum = 10e-12;
tau_bwnum = 1e-12;
gnum = 285; 
alpha_inum = 1.25; 
alpha_mnum = 16; 
ngnum = 3.42;
constnum = 3e10;
format long
Nqwnum = vpa(subs(sol.Nqw,[Isch q V Ab Bb Aqw Bqw Cqw tau_e tau_bw g alpha_i alpha_m ng const],[Ischnum qnum Vnum Abnum Bbnum Aqwnum Bqwnum Cqwnum tau_enum tau_bwnum gnum alpha_inum alpha_mnum ngnum constnum]))

Nqwnum = 

Nbnum = vpa(subs(sol.Nb,[Isch q V Ab Bb Aqw Bqw Cqw tau_e tau_bw g alpha_i alpha_m ng const],[Ischnum qnum Vnum Abnum Bbnum Aqwnum Bqwnum Cqwnum tau_enum tau_bwnum gnum alpha_inum alpha_mnum ngnum constnum]))

Nbnum = 

Npnum = vpa(subs(sol.Np,[Isch q V Ab Bb Aqw Bqw Cqw tau_e tau_bw g alpha_i alpha_m ng const],[Ischnum qnum Vnum Abnum Bbnum Aqwnum Bqwnum Cqwnum tau_enum tau_bwnum gnum alpha_inum alpha_mnum ngnum constnum]))

Npnum = 

for i = 1:6
  error_eqn1(i) = vpa(subs(lhs(eqn1),[Nqw,Nb,Np,Isch q V Ab Bb Aqw Bqw Cqw tau_e tau_bw g alpha_i alpha_m ng const],[Nqwnum(i) Nbnum(i) Npnum(i),Ischnum qnum Vnum Abnum Bbnum Aqwnum Bqwnum Cqwnum tau_enum tau_bwnum gnum alpha_inum alpha_mnum ngnum constnum]));
  error_eqn2(i) = vpa(subs(lhs(eqn2),[Nqw,Nb,Np,Isch q V Ab Bb Aqw Bqw Cqw tau_e tau_bw g alpha_i alpha_m ng const],[Nqwnum(i) Nbnum(i) Npnum(i),Ischnum qnum Vnum Abnum Bbnum Aqwnum Bqwnum Cqwnum tau_enum tau_bwnum gnum alpha_inum alpha_mnum ngnum constnum]));
  error_eqn3(i) = vpa(subs(lhs(eqn3),[Nqw,Nb,Np,Isch q V Ab Bb Aqw Bqw Cqw tau_e tau_bw g alpha_i alpha_m ng const],[Nqwnum(i) Nbnum(i) Npnum(i),Ischnum qnum Vnum Abnum Bbnum Aqwnum Bqwnum Cqwnum tau_enum tau_bwnum gnum alpha_inum alpha_mnum ngnum constnum]));
end
error_eqn1.'

ans = 

error_eqn2.'

ans = 

error_eqn3.'

ans = 

[Mohammad Fozlay]

error_eqn1.' does this show the error for solving the Nqw?

[Torsten]

It shows the value you get if you insert the solutions for Nqw, Nb and Np in the first equation (for the six solutions given).

Answer 2

[Walter Roberson]

syms k [1 3] positive
eqns = Rate_eqns(k)

eqns = 

exact_solution = solve(eqns)

exact_solution = struct with fields:
    k1: (290142196707511*root(z^6 + (7254280273179543825741879705600000000*z^5)/290142196707511 + (125195526563543628214742238301547719652968742878568477720384287106013989665347373648314368*z^4)…
    k2: root(z^6 + (7254280273179543825741879705600000000*z^5)/290142196707511 + (125195526563543628214742238301547719652968742878568477720384287106013989665347373648314368*z^4)/6007463538527559…
    k3: 624173656866248779296875/37778931862957161709568

subs([k1, k2, k3], exact_solution)

ans = 

vpa(ans)

ans = 

function fcns = Rate_eqns(k)

Nqw = k(1);
Nb = k(2);
Np = k(3);

% Ninj = 2e18;
Isch = 15.8e-3; q = 1.6e-19; V = 1e-10;
Ninj = Isch/(q*V);
Ab = 1e8; Bb = 1.2e-10;
tau_b = 1/(Ab+Bb*Nb);
Aqw = 1e8; Bqw = 1.5e-10; Cqw = 5e-30;
tau_r = 1/(Bqw*Nqw); tau_nr = 1/(Aqw+Cqw*Nqw^2); 
tau_e = 10e-12;
tau_bw = 1e-12;
g = 285; alpha_i = 1.25; alpha_m = 16; ng = 3.42;
vg = 3e10/ng;
tau_p = 1/(vg*(alpha_i+alpha_m));


fcns(1) = Ninj - (Nb/tau_b)- (Nb/tau_bw) + (Nqw/tau_e);
fcns(2) = (Nb/tau_bw) -(Nqw/tau_r)- (Nqw/tau_nr)- (Nqw/tau_e) - vg*g*Np;
fcns(3) = vg*g - (Np/tau_p);

end

[Walter Roberson]

There are additional solutions. In three of those solutions, k2 is negative; in the remaining 2 solutions, k2 is complex-valued.

The question does not restrict the range of values of the parameters.

[Mohammad Fozlay]

The Nqw, Nb and Np are all real values as they indicate carrier concentration.

[Walter Roberson]

syms k [1 3] positive
syms R
eqns = Rate_eqns(k)

eqns = 

exact_solution = solve(eqns)

exact_solution = struct with fields:
    k1: (290142196707511*root(z^6 + (7254280273179543825741879705600000000*z^5)/290142196707511 + (125195526563543628214742238301547719652968742878568477720384287106013989665347373648314368*z^4)…
    k2: root(z^6 + (7254280273179543825741879705600000000*z^5)/290142196707511 + (125195526563543628214742238301547719652968742878568477720384287106013989665347373648314368*z^4)/6007463538527559…
    k3: 624173656866248779296875/37778931862957161709568

temp = subs([k1, k2, k3], exact_solution)

temp = 

ur = unique(findSymType(temp, 'root'))

ur = 

wanted_formula = subs(temp, ur, R)

wanted_formula = 

format long g
C = coeffs(children(ur, 1), children(ur,2), 'all')

C = 

DC = roots(double(C))

DC = 
      -8.39273179323287e+21 +  1.02194394498927e+20i
      -8.39273179323287e+21 -  1.02194394498927e+20i
      -8.21401462366212e+21 +                     0i
      -2.92547573165871e+18 +                     0i
      -3.12544333941884e+17 +                     0i
       2.16230193480157e+17 +                     0i

DC = DC(imag(DC) == 0 & real(DC) > 0)

DC = 
      2.16230193480157e+17

vpa(subs(wanted_formula, R, DC))

ans = 

function fcns = Rate_eqns(k)

Nqw = k(1);
Nb = k(2);
Np = k(3);

% Ninj = 2e18;
Isch = 15.8e-3; q = 1.6e-19; V = 1e-10;
Ninj = Isch/(q*V);
Ab = 1e8; Bb = 1.2e-10;
tau_b = 1/(Ab+Bb*Nb);
Aqw = 1e8; Bqw = 1.5e-10; Cqw = 5e-30;
tau_r = 1/(Bqw*Nqw); tau_nr = 1/(Aqw+Cqw*Nqw^2); 
tau_e = 10e-12;
tau_bw = 1e-12;
g = 285; alpha_i = 1.25; alpha_m = 16; ng = 3.42;
vg = 3e10/ng;
tau_p = 1/(vg*(alpha_i+alpha_m));


fcns(1) = Ninj - (Nb/tau_b)- (Nb/tau_bw) + (Nqw/tau_e);
fcns(2) = (Nb/tau_bw) -(Nqw/tau_r)- (Nqw/tau_nr)- (Nqw/tau_e) - vg*g*Np;
fcns(3) = vg*g - (Np/tau_p);

end

What this is telling you is that if you had the list of coefficients shown in C, and the formula in wanted_formula, you could calculate the solutions without the symbolic toolbox.

... but calculating the proper coefficients is probably impractical without the symbolic toolbox.

[Torsten]

@Walter Roberson

Better use

DC = roots(vpa(C))

instead of

DC = roots(double(C))

Otherwise, you'll get a high residual for the second equation:

vpa(subs(eqns,[k(1) k(2) k(3)],vpa(subs(wanted_formula, R, DC))))

[Walter Roberson]

@Torsten

the problem with roots(vpa(C )) is that you get root(expression, z,1), root(expression,z,2) and so on, rather than numeric results. If you are going to vpa() all of that, then you might as well just vpasolve() the original equations.

[Mohammad Fozlay]

M = double(vpa(subs(wanted_formula, R, DC)));

err = double(Rate_eqns(M))

35184000000000	-41304000000000	-0.00048828
%% If I put these solution back into the functions it doesn't give me any converged value

[Torsten]

I don't understand what you are doing there.

The command to get the error in the equations in @Walter Roberson 's code is like I wrote above:

vpa(subs(eqns,[k(1) k(2) k(3)],vpa(subs(wanted_formula, R, DC))))

And are the values for the three unknowns as expected ? :-)

[Mohammad Fozlay]

With your following codes the error seems high for the second variable-

ERR =

[0.000000000037061909097246825695037841796875, -2857595864007.3389535253245412605, 0]

Also the solution when putting those back into the rate equations should yield close to 0 as we formulated our equations as one side is equal to zero. My query is how come the solution is giving large output when we consider them as the inputs for the function?

Thanks

[Torsten]

Read my comment to Walter's code and change

DC = roots(double(C))

to

DC = roots(vpa(C))

Or use my code from above.