how to store the value of intersection points from fzero into a vector

2 vues (au cours des 30 derniers jours)
Dhawal Beohar
Dhawal Beohar le 6 Août 2022
Commenté : Torsten le 7 Août 2022
I am calculating fzero using fzero function. For each value of 'a' (ranging from 0.1 to 1) I need to get value of u*(intersection of fun1 and fun2) as shown in below graph(just for reference). And then i want to store the value of u* in a different vector. Can someone please help me to get the values and store them.Below is my current code which is calculating u* only for a=0.4.
function y = difference(u,d1,n,a,m,T,PsByN_0,UmaxN_0)
d1=20;
n=10^-11.4;
ne=0.5;
m=2.7;
a=0.4;
T=1;
PsByN_0dB=5;
PsByN_0=10.^(PsByN_0dB/10);
UmaxdB = 5;
UmaxN_0=10.^(UmaxdB/10);
u=0.01:0.001:2;
fun1 = @(u) (-1./u).*log(((d1^m)./(a*ne*PsByN_0*T*u+d1^m)*a)./(1-a));
fun2 = @(u) (1./u).*log(((-exp(u.*UmaxN_0).*(exp(-PsByN_0.*u)))./(u.*UmaxN_0+PsByN_0.*u)).*(PsByN_0.*u)-(PsByN_0.*u.*(exp(-PsByN_0.*u))).*(expint(u.*UmaxN_0+PsByN_0.*u))+(exp(-PsByN_0.*u))+((PsByN_0.*u).*(exp(-PsByN_0.*u))).*(expint(PsByN_0.*u))+(exp(u.*UmaxN_0))./((UmaxN_0/PsByN_0)+1));
fun = @(u) (fun1(u) - fun2(u));
%g0=fzero(fun,[0.001,0.02])
g0 = fzero(fun,[0.01])
fun1(g0)
fun2(g0)
plot(u,fun1(u));
hold on;
grid on;
plot(u,fun2(u));
hold on;
grid on;

Réponse acceptée

Torsten
Torsten le 7 Août 2022
Modifié(e) : Torsten le 7 Août 2022
Before using fzero or fsolve, you should check whether fun1 - fun2 really has a zero.
To do this, uncomment the lines
u=[0.01:0.01:2];
plot(u,fun(u,a(3)))
and insert the a-value (in this case a(3)) for which you want to perform the test.
For a>0.5, it seems that fsolve has problems to find a solution. Check it with the method I suggested.
d1=20;
n=10^-11.4;
ne=0.5;
m=2.7;
a=0.01:0.01:0.5;
T=1;
PsByN_0dB=5;
PsByN_0=10.^(PsByN_0dB/10);
UmaxdB = 5;
UmaxN_0=10.^(UmaxdB/10);
fun1 = @(u,a) (-1./u).*log(((d1.^m)./(a.*ne.*PsByN_0.*T.*u+d1.^m).*a)./(1-a));
fun2 = @(u) (1./u).*log(((-exp(u.*UmaxN_0).*(exp(-PsByN_0.*u)))./(u.*UmaxN_0+PsByN_0.*u)).*(PsByN_0.*u)-(PsByN_0.*u.*(exp(-PsByN_0.*u))).*(expint(u.*UmaxN_0+PsByN_0.*u))+(exp(-PsByN_0.*u))+((PsByN_0.*u).*(exp(-PsByN_0.*u))).*(expint(PsByN_0.*u))+(exp(u.*UmaxN_0))./((UmaxN_0/PsByN_0)+1));
fun = @(u,a) (fun1(u,a) - fun2(u));
%u=[0.01:0.01:2];
%plot(u,fun(u,a(3)))
%g0=fzero(fun,[0.001,0.02])
options = optimset('Display','none');
g0 = arrayfun(@(a)fsolve(@(u)fun(u,a),[0.01],options),a)
g0 = 1×50
1.6739 1.4530 1.3231 1.2304 1.1581 1.0986 1.0480 1.0039 0.9646 0.9293 0.8971 0.8674 0.8399 0.8141 0.7900 0.7671 0.7455 0.7248 0.7050 0.6861 0.6678 0.6501 0.6330 0.6164 0.6001 0.5843 0.5688 0.5536 0.5386 0.5238
%fun(g0,a)
plot(a,g0)
  5 commentaires
Torsten
Torsten le 7 Août 2022
Modifié(e) : Torsten le 7 Août 2022
Alright....So i don't need to use any for loop or something to iterate through the values of a and then store the value of u* in other vector?
This is all done in the line
g0 = arrayfun(@(a)fsolve(@(u)fun(u,a),[0.01],options),a)
If you have more experience with for-loops: this code is equivalent to the one above, but you have more flexibility providing the initial guess for u*:
d1=20;
n=10^(-11.4);
ne=0.5;
m=2.7;
a=0.01:0.01:0.5;
T=1;
PsByN_0dB=5;
PsByN_0=10.^(PsByN_0dB/10);
UmaxdB = 5;
UmaxN_0=10.^(UmaxdB/10);
fun1 = @(u,a) (-1./u).*log(((d1.^m)./(a.*ne.*PsByN_0.*T.*u+d1.^m).*a)./(1-a));
fun2 = @(u) (1./u).*log(((-exp(u.*UmaxN_0).*(exp(-PsByN_0.*u)))./(u.*UmaxN_0+PsByN_0.*u)).*(PsByN_0.*u)-(PsByN_0.*u.*(exp(-PsByN_0.*u))).*(expint(u.*UmaxN_0+PsByN_0.*u))+(exp(-PsByN_0.*u))+((PsByN_0.*u).*(exp(-PsByN_0.*u))).*(expint(PsByN_0.*u))+(exp(u.*UmaxN_0))./((UmaxN_0/PsByN_0)+1));
fun = @(u,a) (fun1(u,a) - fun2(u));
options = optimset('Display','none');
g00 = 0.01;
g0 = zeros(size(a));
for i=1:numel(a)
g0(i) = fsolve(@(u)fun(u,a(i)),g00,options);
g00 = g0(i);
end
plot(a,g0)
Torsten
Torsten le 7 Août 2022
@Dhawal Beohar comment moved here
Thanks for explaining things. I really appreciate your help. Thank you so much for being there !!!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Optimization dans Help Center et File Exchange

Produits


Version

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by