Different results from directly using fsolve function and using the code within program in fsolve setting

3 visualizaciones (últimos 30 días)
i want to use the fsolve commad to solve the roots of four variables, during the calcuation, the return flag always be -2.
so i change manys times about the start point as shown below.
in one partical setting, the model genetate value(na_fa_medel na_fr_model etc..) close to the target (na_fa na_fr etc..) but the flag still be -2.
so i seperately run the fslove program and the program within the fsolve function
the return value from the former situation is [0.0021997 0.012094 0.14924 0.015243] , howver the ture value should close to [0.0022 0.0121 0.0171 0.0087] based on later situation.
in these four values, it seems that the former two items is correct but the latter two is incorrect.
i don't know why there is such a difference but i have checked many times that the code is correct.
can anyone helps me to fix it or tell me why?
if you directly use the fsolve founction and fval, the results is [0.0021997 0.012094 0.14924 0.015243], however , if you run the program within the program , the value should close to [0.0022 0.0121 0.0171 0.0087]
Here is the program:
clear all;
close all;
clc;
format short g
global Na Nr Nu Na_fa Na_pr Na_pu Na_fr Na_fu Aa Ar wr Pa
Na_fa =0.206;
Na_pr = 0.197;
Na_pu = 0.29;
Na_fr = 0.127;
Na_fu = 0.18;
Ar = 1;%normalization
wr = Ar;%defination
MU_vec = [0 0 0 0];
SIGMA_a = 0.4646;
SIGMA_phi = 1.1687;
SIGMA_r = 0.8127;
SIGMA_u = 0.4394;
COV_ap = -0.4713;
COV_ar = 0.0154;
COV_au = 0.0052;
COV_pr = 0; %based on defination
COV_pu = 0; %based on defination
COV_ru = 0; %based on defination
SIGMA_mat = [SIGMA_a^2 COV_ap COV_ar COV_au; COV_ap SIGMA_phi^2 COV_pr COV_pu; COV_ar COV_pr SIGMA_r^2 COV_ru; COV_au COV_pu COV_ru SIGMA_u^2];
n = 1000000;
Ngen = mvnrnd(MU_vec,SIGMA_mat,n);
LNgen = exp(Ngen);
vecSA = LNgen(:,1);
vecFI = LNgen(:,2);
vecSR = LNgen(:,3);
vecSU = LNgen(:,4);
ETA_R = 0.6019;
ETA_U = 0.6641;
wa = 0.3981;
wu = 1.1863;
options = optimset('Display','iter','TolX',1e-10,'LargeScale','on');
cr = 0.075;%assume
cu = 0.1;%assume
kappaa = 1.25;%assume
kappau = 1.3;%assume
nuv =0.4;
guess = [kappaa, kappau, cr, cu];
[x, fval, exitflag,~] = fsolve(@solve_kappaandc,guess,options,vecSA,vecFI,vecSR,vecSU,nuv,wa,wu,ETA_R,ETA_U);
fval;
abs(fval);
The program is:;
function f = solve_kappaandc(x,BARRIER_R,BARRIER_U,vecSA,vecFI,vecSR,vecSU,nuv,wa,wu,ETA_R,ETA_U);
global wr Na_fr Na_fu Na_pr Na_fa n Na_pu ;
Na_fa =0.206;
Na_pr = 0.197;
Na_pu = 0.29;
Na_fr = 0.127;
Na_fu = 0.18;
kappaa = x(1);
kappau = x(2);
cr = x(3);
cu = x(4);
Income_a = kappaa.*wa.*vecSA.*vecFI;
Income_fr = (1-ETA_R)*wr.*vecSR;
Income_fu = (1-ETA_U)*kappau.*wu.*vecSU;
vce_cr = ones(n,1)*cr;
vce_cu = ones(n,1)*cu;
oneminuscr = 1-vce_cr;
oneminuscu = 1-vce_cu;
time_pr = (Income_a*nuv./Income_fr).^(1/(1-nuv));
selec = [time_pr oneminuscr];
time_pr = min(selec,[],2);
time_pu = (Income_a*nuv./Income_fu).^(1/(1-nuv));
selec = [time_pu oneminuscu];
time_pu = min(selec,[],2);
Income_pr = (1-cr-time_pr).*Income_fr + Income_a.*time_pr.^nuv;
Income_pu = (1-cu-time_pu).*Income_fu + Income_a.*time_pu.^nuv;
OC_a = zeros(n,1);
OC_a(Income_a > Income_fr & Income_a>Income_fu & Income_a>Income_pr & Income_a>Income_pu)=1;
OC_fr = zeros(n,1);
OC_fr(Income_fr > Income_a & Income_fr>Income_fu & Income_fr>Income_pr & Income_fr>Income_pu)=1;
OC_fu = zeros(n,1);
OC_fu(Income_fu > Income_a & Income_fu>Income_fr & Income_fu>Income_pr & Income_fu>Income_pu)=1;
OC_pr = zeros(n,1);
OC_pr(Income_pr > Income_a & Income_pr>Income_fr & Income_pr>Income_fu & Income_pr>Income_pu)=1;
OC_pu = zeros(n,1);
OC_pu(Income_pu > Income_a & Income_pu>Income_fr & Income_pu>Income_fu & Income_pu>Income_pr)=1;
N = sum(ones(length(OC_a),1));
Na_fa_model = sum(OC_a)/N;
Na_fr_model = sum(OC_fr)/N;
Na_fu_model = sum(OC_fu)/N;
Na_pr_model = sum(OC_pr)/N;
Na_pu_model = sum(OC_pu)/N;
f1 = Na_fr_model - Na_fr;
f2 = Na_fu_model - Na_fu;
f3 = Na_pr_model - Na_pr;
f4 = Na_pu_model - Na_pu;
f = [f1 f2 f3 f4 ];

Respuesta aceptada

Dyuman Joshi
Dyuman Joshi el 2 de Jun. de 2023
Using global is generally not recommended, specially when it is not required.
There are many errors in your code, I have edited them.
format short g
Ar = 1;%normalization
wr = Ar;%defination
MU_vec = [0 0 0 0];
SIGMA_a = 0.4646;
SIGMA_phi = 1.1687;
SIGMA_r = 0.8127;
SIGMA_u = 0.4394;
COV_ap = -0.4713;
COV_ar = 0.0154;
COV_au = 0.0052;
COV_pr = 0; %based on defination
COV_pu = 0; %based on defination
COV_ru = 0; %based on defination
SIGMA_mat = [SIGMA_a^2 COV_ap COV_ar COV_au; COV_ap SIGMA_phi^2 COV_pr COV_pu; COV_ar COV_pr SIGMA_r^2 COV_ru; COV_au COV_pu COV_ru SIGMA_u^2];
n = 1000000;
Ngen = mvnrnd(MU_vec,SIGMA_mat,n);
LNgen = exp(Ngen);
vecSA = LNgen(:,1);
vecFI = LNgen(:,2);
vecSR = LNgen(:,3);
vecSU = LNgen(:,4);
ETA_R = 0.6019;
ETA_U = 0.6641;
wa = 0.3981;
wu = 1.1863;
options = optimset('Display','iter','TolX',1e-10,'LargeScale','on');
cr = 0.075;%assume
cu = 0.1;%assume
kappaa = 1.25;%assume
kappau = 1.3;%assume
nuv =0.4;
guess = [kappaa, kappau, cr, cu];
%You can either define all the variables inside the function to be solved
%Or you can pass them as input arguements (which is what I have done here)
[x, fval, exitflag] = fsolve(@(x) solve_kappaandc(x,vecSA,vecFI,vecSR,vecSU,nuv,wa,wu,ETA_R,ETA_U,wr,n),guess,options)
Norm of First-order Trust-region Iteration Func-count ||f(x)||^2 step optimality radius 0 5 0.000503699 0 1 Equation solved at initial point. fsolve completed because the vector of function values at the initial point is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
x = 1×4
1.25 1.3 0.075 0.1
fval = 1×4
0.001515 -0.011287 0.01714 0.008957
exitflag =
1
function f = solve_kappaandc(x,vecSA,vecFI,vecSR,vecSU,nuv,wa,wu,ETA_R,ETA_U,wr,n);
Na_fa =0.206;
Na_pr = 0.197;
Na_pu = 0.29;
Na_fr = 0.127;
Na_fu = 0.18;
kappaa = x(1);
kappau = x(2);
cr = x(3);
cu = x(4);
Income_a = kappaa.*wa.*vecSA.*vecFI;
Income_fr = (1-ETA_R)*wr.*vecSR;
Income_fu = (1-ETA_U)*kappau.*wu.*vecSU;
vce_cr = ones(n,1)*cr;
vce_cu = ones(n,1)*cu;
oneminuscr = 1-vce_cr;
oneminuscu = 1-vce_cu;
time_pr = (Income_a*nuv./Income_fr).^(1/(1-nuv));
selec = [time_pr oneminuscr];
time_pr = min(selec,[],2);
time_pu = (Income_a*nuv./Income_fu).^(1/(1-nuv));
selec = [time_pu oneminuscu];
time_pu = min(selec,[],2);
Income_pr = (1-cr-time_pr).*Income_fr + Income_a.*time_pr.^nuv;
Income_pu = (1-cu-time_pu).*Income_fu + Income_a.*time_pu.^nuv;
OC_a = zeros(n,1);
OC_a(Income_a > Income_fr & Income_a>Income_fu & Income_a>Income_pr & Income_a>Income_pu)=1;
OC_fr = zeros(n,1);
OC_fr(Income_fr > Income_a & Income_fr>Income_fu & Income_fr>Income_pr & Income_fr>Income_pu)=1;
OC_fu = zeros(n,1);
OC_fu(Income_fu > Income_a & Income_fu>Income_fr & Income_fu>Income_pr & Income_fu>Income_pu)=1;
OC_pr = zeros(n,1);
OC_pr(Income_pr > Income_a & Income_pr>Income_fr & Income_pr>Income_fu & Income_pr>Income_pu)=1;
OC_pu = zeros(n,1);
OC_pu(Income_pu > Income_a & Income_pu>Income_fr & Income_pu>Income_fu & Income_pu>Income_pr)=1;
N = sum(ones(length(OC_a),1));
Na_fa_model = sum(OC_a)/N;
Na_fr_model = sum(OC_fr)/N;
Na_fu_model = sum(OC_fu)/N;
Na_pr_model = sum(OC_pr)/N;
Na_pu_model = sum(OC_pu)/N;
f1 = Na_fr_model - Na_fr;
f2 = Na_fu_model - Na_fu;
f3 = Na_pr_model - Na_pr;
f4 = Na_pu_model - Na_pu;
f = [f1 f2 f3 f4 ];
end
  2 comentarios
Peinan Peinan
Peinan Peinan el 2 de Jun. de 2023
Thanks a lot, i have benn tormented by this problem for almost two days and i basiclly not think about that the problem comes from the global.
Really thank you !!!!!!!
Learned a lot from your answer!!
Dyuman Joshi
Dyuman Joshi el 2 de Jun. de 2023
There were errors as well other than the use of global, notably how to define the function handle as input to fsolve()
Compare the lines where fsolve() is called in initial code and the corrected code to observe the difference.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by