'fsolve' producing error message
显示 更早的评论
Hi the community,
I've been struggling with the following code for a few days now, not sure what is wrong with it. I would appreciate if you have a look at the code. It is very simple code, a simulation exercise. I am getting the following message
- Error using trustnleqn (line 28)
- Objective function is returning undefined values at initial point. FSOLVE cannot continue.
- Error in fsolve (line 366)
- [x,FVAL,JACOB,EXITFLAG,OUTPUT,msgData]=...
- Error in paper2_v61 (line 76)
- out1=fsolve(@productivity,x0,options);
Thank you!
The code: clc clear all
global A tauc tauw a b sigma alpha S phic phiw rc rw
format short
T=200;
alpha=0.75;
tauw=1;
tauc=2;
rc=3;
rw=6;
sigma=3;
phic=0.5;
phiw=1;
a=1;
b=0.81;
%Initial values
Ac0=1;
Aw0=5;
Sw0=7000;
Sc0=35000;
Scb=10^6;
Theta=zeros(T,1);
Ec=zeros(T,1);
Ew=zeros(T,1);
X=zeros(T,1);
Yc=zeros(T,1);
Yw=zeros(T,1);
Ac=zeros(T+1,1);
Aw=zeros(T+1,1);
Ac(1,1)=Ac0;
Aw(1,1)=Aw0;
Pc=zeros(T,1);
Pw=zeros(T,1);
pc=zeros(T,1);
pw=zeros(T,1);
cc=zeros(T,1);
cw=zeros(T,1);
Nc=zeros(T,1);
Nw=zeros(T,1);
Sc=zeros(T+1,1);
Sw=zeros(T+1,1);
Sc(1,1)=Sc0;
Sw(1,1)=Sw0;
x0=[0.5,0.5,100];
for t=1:T
A=[Ac(t,1),Aw(t,1)];
S=[Sc(t,1),Sw(t,1)];
options=optimset('TolX',1e-9,'TolFun',1e-9);
out1=fsolve(@productivity,x0,options);
Nc(t,1)=out1(2);
Nw(t,1)=1-out1(2);
Theta(t,1)=out1(1);
Pc(t,1)=out1(3);
Ac(t+1,1)=(1+tauc*Theta(t,1))*Ac(t,1);
Aw(t+1,1)=(1+tauw*Theta(t,1))*Aw(t,1);
X(t,1)=Sc(t,1)*(1/Ac(t+1,1));
Sc(t,1)=Sc(t,1)+X(t,1);
Pw(t,1)=(1-Pc(t,1)^(1-sigma))^(1/(1-sigma));
cc(t,1)=Pc(t,1)^(1/(1-alpha))*Ac(t+1,1)*(1-alpha);
cw(t,1)=Pw(t,1)^(1/(1-alpha))*Aw(t+1,1)*(1-alpha);
pc(t,1)=cc(t,1);
pw(t,1)=cw(t,1);
Ec(t,1)=(1-(cc(t,1)/phic)^1/rc)*Sc(t,1);
Ew(t,1)=(1-(cw(t,1)/phiw)^1/rw)*Sw(t,1);
Sc(t+1)=Sc(t,1)-Ec(t,1);
Sw(t+1)=Sw(t,1)-Ew(t,1);
end
Function
function F=productivity(x)
global a b tauc tauw A alpha S sigma phiw phic rc rw
Ac=(1+tauc*x(1))*A(1);
Aw=(1+tauc*x(1))*A(2);
Sc=S(1);
Sw=S(2);
cc=x(3)^(1/(1-alpha))*Ac*(1-alpha);
cw=((1-x(3)^(1-sigma))^(1/(1-sigma)))^(1/(1-alpha))*Aw*(1-alpha);
E1=((1-(cc/phic)^1/rc)*Sc)/((1-(cw/phiw)^1/rw)*Sw);
E2=(x(3)/((1-x(3)^(1-sigma))^(1/(1-sigma))))^(alpha/(alpha-1))*(Ac/Aw)^(sigma*(1-alpha)-1)*(cc/cw)^-sigma;
F=[ -x(1)+a*(x(2)/((1+tauc)*A(1)))^b;
-x(1)+a*((1-x(2))/((1+tauw)*A(2)))^b;
E1-E2;
];
采纳的回答
Alan Weiss
2014-12-16
You didn't include your call to fsolve so it is hard to know exactly what you are doing. We don't know your options, your initial point x0, or your fsolve call.
Nevertheless, the error is clear. You have an initial point for fsolve, I suppose it is x0. When evaluating
productivity(x0)
fsolve found that there were NaN in the result. Either choose an initial point that does not produce this behavior, or debug your program so that it does not produce NaN results.
Alan Weiss
MATLAB mathematical toolbox documentation
更多回答（1 个）
When the loop reaches t=46, I find the following
K>> productivity(x0)
ans =
-0.0917
-0.4998
Inf
I assume the problem is now clear. productivity() is not returning finite real values at the initial point x0. It is because Sc has reached a super-huge number
K>> Sc
Sc =
3.2735e+301
making E1 numerically infinite
K>> E1
E1 =
Inf
类别
在 帮助中心 和 File Exchange 中查找有关 Digital Filter Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
