Different results from directly using fsolve function and using the code within program in fsolve setting
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Peinan Peinan
el 2 de Jun. de 2023
Comentada: Dyuman Joshi
el 2 de Jun. de 2023
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 ];
0 comentarios
Respuesta aceptada
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)
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
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.
Más respuestas (0)
Ver también
Categorías
Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!