Matrix input in fzero solver

4 次查看(过去 30 天)
Hello,
I am stuck with my code. I am trying to vectorize loops for performance issues. However, I have a problem when the vectorization involves fzero. Here I post an example of what I want to get. First code line shows the script with the loops, and the second one the function file.
MAIN SCRIPT
clear all
clc
grida=10;
gridz=60;
w=1.296;
mu1=0.85;
theta=0.1;
s = rng;
bbar_AQ = zeros(grida,gridz);
pistar_n_f = zeros(1,gridz);
k_star_f = zeros(1,gridz);
auxiliar = zeros(grida,gridz);
debt = zeros(grida,gridz);
A=randn(10,60);
aa=0.001:(w-0.002)/(grida-1):w-0.001;
z = linspace(1, 12, gridz);
opt=optimset('display', 'off');
for j=1:grida
for t=1:gridz
k_star_f(t) =(z(t))^(1/(1-mu1-theta));
pistar_n_f(t) = (1-mu1)*(mu1/w)^(mu1/(1-mu1))*(z(t))^(1/(1-mu1))*(k_star_f(t)^(theta/(1-mu1)));
bbar_AQ(j,t) = (pistar_n_f(t) + aa(j))/0.75;
auxiliar(j,t) =k_star_f(t);
if bbar_AQ(j,t)>=auxiliar(j,t)
debt(j,t)=A(j,t);
else
if j==1
bbar_0 = 0.1;
else
bbar_0 = 0.11;
end
debt(j,t)=fzero(@(x)b_max(x,z(t),aa(j),theta,mu1,w),bbar_0,opt);
end
end
end
FUNCTION FILE
function f=b_max(x,z,a,theta,mu,w)
f=(((1-mu)*(mu/w)^(mu/(1-mu))*(z*(x+a)^theta)^(1/(1-mu))))-x;
The above code runs fine. Howeve, I am trying to vectorize it but when I arrive to the fzero line of code, I am stuck about how to proceed. This is what I did:
clear all
clc
grida = 10;
gridz = 60;
w = 1.296;
mu1 = 0.85;
theta = 0.1;
aa = 0.001:(w-0.002)/(grida-1):w-0.001;
z = linspace(1, 12, gridz);
rng(s)
A=randn(10,60);
% Compute k_star_f and pistar_n_f using vectorized operations
k_star_f = z.^(1./(1-mu1-theta));
pistar_n_f = (1-mu1)*(mu1/w).^(mu1./(1-mu1)).*(z.^(1./(1-mu1))).*(k_star_f.^(theta./(1-mu1)));
% Compute bbar_AQ using vectorized operations
expanded_aa = repmat(aa(:), 1, gridz);
bbar_AQ = (pistar_n_f + expanded_aa) / 0.75;
% Compute auxiliar using vectorized operations
auxiliar = k_star_f;
auxiliar=repmat(auxiliar,grida,1);
% Compute debt using vectorized operations based on the condition
debt2 = A2 .* (bbar_AQ2 >= auxiliar2) + fzero(@(x)b_max(x,z(t),aa(j),theta,mu1,w),bbar_0,opt).* (bbar_AQ2 < auxiliar2);
Any help would be really appreciated.
Thank in advanced!
  8 个评论
Torsten
Torsten 2023-9-8
编辑:Torsten 2023-9-8
So, it will not possible to use parfor loop?
No. The loop as written has to be processed sequentially because you made the initial values of the i-th call to "fzero" depend on the result of the (i-1)th call.
You might try to set bbar_0 outside and only leave the line
bbar_constrained(j,t)=fzero(@(x)b_max_2_OnlyPi(x,z(t),eta,aa(j),theta,mu1,w),bbar_0,opt);
within the double loop.
This loop should then be possible to be processed in parallel.
Ibai Ostolozaga Falcon
Thank you so much for your help Torsten!

请先登录,再进行评论。

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Loops and Conditional Statements 的更多信息

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by