从 Solve 到 Fsolve
From Solve to Fsolve
我想解三个方程式,三个未知数。
我用 symbolic toolbox
指定方程式。我知道我可以使用 solve
函数让 matlab 帮我找到一个数值解。但是,对于 3 个未知数中的 3 个方程,matlab 应该能够找到解析解 (fsolve
)。我只是不确定如何更改代码以便我可以使用 fsolve
而不是 solve
.
在我的代码下面:
全部清除
syms Kl Kh alpha nu w phi delta P beta zh zl Ezh Ezl
nu1 = (1/(1-nu));
f1 = ((zl * (Kl^alpha))^nu1 + (zh * (Kh^alpha))^nu1) * nu^(nu*nu1) * (w^(-nu*nu1)) - w/phi + delta*(Kl + Kh)*P
f2 = Kh - (( (1-beta*(1-delta))*P * (w^(nu1*nu))*(nu^(nu*nu1)) ) / (beta*alpha* (Ezh)^nu1) )^((1-nu)/(alpha+nu-1))
f3 = Kl - (( (1-beta*(1-delta))*P * (w^(nu1*nu))*(nu^(nu*nu1)) ) / (beta*alpha* (Ezl)^nu1) )^((1-nu)/(alpha+nu-1))
f1 = subs(f1, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95})
f2 = subs(f2, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95})
f3 = subs(f3, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95})
S = solve([f1 == 0, f2 == 0, f3 == 0],...
[w, Kh, Kl], 'ReturnConditions', true);
同时我找到了解决方案。
这是完成代码:
function SSfunction = SSfunction(x)
syms alphaa nu phi delta p betaa zh zl ezh ezl
nu1 = (1/(1-nu));
f1 = ((zl * (x(3)^alphaa))^nu1 + (zh * (x(2)^alphaa))^nu1) * nu^(nu*nu1) * (x(1)^(-nu*nu1)) - x(1)/phi - delta*(x(3) + x(2))*p;
f2 = x(2) - (( betaa*alphaa*(ezh^(nu1)) * (nu^(nu*nu1)) )/( (1-betaa*(1-delta))*p* (x(1)^(nu*nu1)) ))^((1-nu)/(1-alphaa-nu));
f3 = x(3) - (( betaa*alphaa*(ezl^(nu1)) * (nu^(nu*nu1)) )/( (1-betaa*(1-delta))*p* (x(1)^(nu*nu1)) ))^((1-nu)/(1-alphaa-nu));
f1 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950});
f3 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950});
f2 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950});
SSfunction(1) = eval(f1)
SSfunction(2) = eval(f2)
SSfunction(3) = eval(f3)
end
x0 = [1,2,0.7];
fun = @SSfunction;
x = fsolve(fun,x0)
使用matlabFunction
将您的符号表达式转换为向量化数值函数,可以直接被fsolve
使用:
...
f = matlabFunction([f1;f2;f3],'Vars',{[w;Kh;Kl]});
w0 = 1;
Kh0 = 1;
Kl0 = 1;
x0 = [w0;Kh0;Kl0];
x = fsolve(f,x0)
这将比使用 fsolve
中的符号表达式本身快一个数量级。为了提高速度,您还可以通过手动矢量化函数来完全摆脱符号数学:
alpha = 0.27;
beta = 0.96;
nu = 0.6;
phi = 2.15;
delta = 0.065;
zh = 1.11687642219068;
zl = 0.895354204038589;
Ezh = 1.07811003137331;
Ezl = 0.934120594855956;
P = 0.95;
nu1 = (1/(1-nu));
f = @(w,Kh,Kl)[((zl * (Kl.^alpha))^nu1 + (zh * (Kh.^alpha))^nu1) * nu^(nu*nu1) .* (w.^(-nu*nu1)) - w/phi + delta*(Kl + Kh)*P;
Kh - (( (1-beta*(1-delta))*P * (w.^(nu1*nu))*(nu^(nu*nu1)) ) / (beta*alpha* (Ezh)^nu1) )^((1-nu)/(alpha+nu-1));
Kl - (( (1-beta*(1-delta))*P * (w.^(nu1*nu))*(nu^(nu*nu1)) ) / (beta*alpha* (Ezl)^nu1) )^((1-nu)/(alpha+nu-1))];
w0 = 1;
Kh0 = 1;
Kl0 = 1;
x0 = [w0;Kh0;Kl0];
x = fsolve(@(x)f(x(1,:),x(2,:),x(3,:)),x0)
我想解三个方程式,三个未知数。
我用 symbolic toolbox
指定方程式。我知道我可以使用 solve
函数让 matlab 帮我找到一个数值解。但是,对于 3 个未知数中的 3 个方程,matlab 应该能够找到解析解 (fsolve
)。我只是不确定如何更改代码以便我可以使用 fsolve
而不是 solve
.
在我的代码下面:
全部清除
syms Kl Kh alpha nu w phi delta P beta zh zl Ezh Ezl
nu1 = (1/(1-nu));
f1 = ((zl * (Kl^alpha))^nu1 + (zh * (Kh^alpha))^nu1) * nu^(nu*nu1) * (w^(-nu*nu1)) - w/phi + delta*(Kl + Kh)*P
f2 = Kh - (( (1-beta*(1-delta))*P * (w^(nu1*nu))*(nu^(nu*nu1)) ) / (beta*alpha* (Ezh)^nu1) )^((1-nu)/(alpha+nu-1))
f3 = Kl - (( (1-beta*(1-delta))*P * (w^(nu1*nu))*(nu^(nu*nu1)) ) / (beta*alpha* (Ezl)^nu1) )^((1-nu)/(alpha+nu-1))
f1 = subs(f1, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95})
f2 = subs(f2, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95})
f3 = subs(f3, {alpha, beta, nu, phi,delta, zh, zl, Ezh, Ezl, P}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.95})
S = solve([f1 == 0, f2 == 0, f3 == 0],...
[w, Kh, Kl], 'ReturnConditions', true);
同时我找到了解决方案。
这是完成代码:
function SSfunction = SSfunction(x)
syms alphaa nu phi delta p betaa zh zl ezh ezl
nu1 = (1/(1-nu));
f1 = ((zl * (x(3)^alphaa))^nu1 + (zh * (x(2)^alphaa))^nu1) * nu^(nu*nu1) * (x(1)^(-nu*nu1)) - x(1)/phi - delta*(x(3) + x(2))*p;
f2 = x(2) - (( betaa*alphaa*(ezh^(nu1)) * (nu^(nu*nu1)) )/( (1-betaa*(1-delta))*p* (x(1)^(nu*nu1)) ))^((1-nu)/(1-alphaa-nu));
f3 = x(3) - (( betaa*alphaa*(ezl^(nu1)) * (nu^(nu*nu1)) )/( (1-betaa*(1-delta))*p* (x(1)^(nu*nu1)) ))^((1-nu)/(1-alphaa-nu));
f1 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950});
f3 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950});
f2 = subs(f1, {alphaa, betaa, nu, phi,delta, zh, zl, ezh, ezl, p}, {0.27, 0.96, 0.60, 2.15,0.065,1.11687642219068,0.895354204038589,1.07811003137331,0.934120594855956, 0.950});
SSfunction(1) = eval(f1)
SSfunction(2) = eval(f2)
SSfunction(3) = eval(f3)
end
x0 = [1,2,0.7];
fun = @SSfunction;
x = fsolve(fun,x0)
使用matlabFunction
将您的符号表达式转换为向量化数值函数,可以直接被fsolve
使用:
...
f = matlabFunction([f1;f2;f3],'Vars',{[w;Kh;Kl]});
w0 = 1;
Kh0 = 1;
Kl0 = 1;
x0 = [w0;Kh0;Kl0];
x = fsolve(f,x0)
这将比使用 fsolve
中的符号表达式本身快一个数量级。为了提高速度,您还可以通过手动矢量化函数来完全摆脱符号数学:
alpha = 0.27;
beta = 0.96;
nu = 0.6;
phi = 2.15;
delta = 0.065;
zh = 1.11687642219068;
zl = 0.895354204038589;
Ezh = 1.07811003137331;
Ezl = 0.934120594855956;
P = 0.95;
nu1 = (1/(1-nu));
f = @(w,Kh,Kl)[((zl * (Kl.^alpha))^nu1 + (zh * (Kh.^alpha))^nu1) * nu^(nu*nu1) .* (w.^(-nu*nu1)) - w/phi + delta*(Kl + Kh)*P;
Kh - (( (1-beta*(1-delta))*P * (w.^(nu1*nu))*(nu^(nu*nu1)) ) / (beta*alpha* (Ezh)^nu1) )^((1-nu)/(alpha+nu-1));
Kl - (( (1-beta*(1-delta))*P * (w.^(nu1*nu))*(nu^(nu*nu1)) ) / (beta*alpha* (Ezl)^nu1) )^((1-nu)/(alpha+nu-1))];
w0 = 1;
Kh0 = 1;
Kl0 = 1;
x0 = [w0;Kh0;Kl0];
x = fsolve(@(x)f(x(1,:),x(2,:),x(3,:)),x0)