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.
2 views (last 30 days)
Show older comments
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);
% Check the exit flag
if exitflag <= 0
error('fsolve did not converge to a solution.')
end
% 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
3 Comments
Torsten
on 25 Apr 2023
Edited: Torsten
on 25 Apr 2023
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.
Accepted Answer
Torsten
on 25 Apr 2023
Edited: Torsten
on 25 Apr 2023
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])
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]))
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]))
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]))
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.'
error_eqn2.'
error_eqn3.'
2 Comments
More Answers (1)
Walter Roberson
on 25 Apr 2023
Moved: Walter Roberson
on 25 Apr 2023
syms k [1 3] positive
eqns = Rate_eqns(k)
exact_solution = solve(eqns)
subs([k1, k2, k3], exact_solution)
vpa(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
9 Comments
Torsten
on 26 Apr 2023
Read my comment to Walter's code and change
DC = roots(double(C))
to
DC = roots(vpa(C))
Or use my code from above.
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!