Is this a correct way to use fsolve?

Question

Salvatore Bianco 2024-7-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2137318-is-this-a-correct-way-to-use-fsolve

编辑： Torsten 2024-7-15

Since fsolve keeps giving me answers with a very small but non-zero imaginary part, which i really don't want, I though about giving the derivative of my function. Is this a correct way to do it? Also, is there a way to tell the function not to go outside the reale line?

%g is a function that is defined by an equation. I know that is invertible
%and takes values between 0 and kappa.
function g = g(y, z, p, kappa, beta, mu, muz, sigma, sigmaz)
K = (mu^2/sigma^2)*0.5;
M = (sigmaz*mu)/sigma;
q = 1/(p-1);
alpha = (sqrt((beta-M-K)^2+4*K*(beta-muz))-beta+M+K)/(2*K);
C = beta-K*p/(1-p);
%function whose zero i need to find (with respect to w)
F = @(w) ((1-p)/C)*(kappa^q - w^q)+y+z-z*(w/kappa)^(alpha-1);
%its derivative with respect to w
J = @(w) ((1-p)/C)*(-q*(w^(q-1))) -z*((alpha-1)*w^(alpha-2))/(kappa^(alpha-1));
options = optimoptions('fsolve', 'SpecifyObjectiveGradient', true);
fun = {F, J};
w0 = kappa/2;
g = fsolve(fun, w0, options);
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-7-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2137318-is-this-a-correct-way-to-use-fsolve#answer_1485408

编辑：Torsten 2024-7-15

在 MATLAB Online 中打开

Solve in w^2 instead of w - then there shouldn't be imaginary parts in the solution:

F = @(w) ((1-p)/C)*(kappa^q - (w^2)^q)+y+z-z*(w^2/kappa)^(alpha-1);

instead of

F = @(w) ((1-p)/C)*(kappa^q - w^q)+y+z-z*(w/kappa)^(alpha-1);

Don't forget to take the square of g before exiting the function:

...
g = fsolve(fun, w0, options);
g = g^2;
end

8 个评论
显示 6更早的评论隐藏 6更早的评论

Torsten 2024-7-15

编辑：Torsten 2024-7-15

在 MATLAB Online 中打开

clear all
n = 100;
y = linspace(0, 10, n);
z = [0, 2, 5, 10];
m = length(z);
beta = 8; kappa = 4;
mu = 1; sigma = 1;
muz = 1; sigmaz = 1;
p = -2;
for i=1:m
for j=1:n
    [G(i, j),error(i,j)]= g(y(j), z(i), p, kappa, beta, mu, muz, sigma, sigmaz);
end
end
G
G = 4x100
    3.9999    1.3246    0.5917    0.3137    0.1859    0.1190    0.0808    0.0573    0.0421    0.0318    0.0247    0.0195    0.0157    0.0128    0.0106    0.0088    0.0075    0.0064    0.0055    0.0047    0.0041    0.0036    0.0032    0.0028    0.0025    0.0022    0.0020    0.0018    0.0016    0.0015
    3.9999    1.3246    0.5917    0.3137    0.1859    0.1190    0.0808    0.0573    0.0421    0.0318    0.0247    0.0195    0.0157    0.0128    0.0106    0.0088    0.0075    0.0064    0.0055    0.0047    0.0041    0.0036    0.0032    0.0028    0.0025    0.0022    0.0020    0.0018    0.0016    0.0015
    3.9999    1.3246    0.5917    0.3137    0.1859    0.1190    0.0808    0.0573    0.0421    0.0318    0.0247    0.0195    0.0157    0.0128    0.0106    0.0088    0.0075    0.0064    0.0055    0.0047    0.0041    0.0036    0.0032    0.0028    0.0025    0.0022    0.0020    0.0018    0.0016    0.0015
    3.9999    1.3246    0.5917    0.3137    0.1859    0.1190    0.0808    0.0573    0.0421    0.0318    0.0247    0.0195    0.0157    0.0128    0.0106    0.0088    0.0075    0.0064    0.0055    0.0047    0.0041    0.0036    0.0032    0.0028    0.0025    0.0022    0.0020    0.0018    0.0016    0.0015
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
max(abs(error(:)))
ans = 2.6811e-06
function [g,error] = g(y, z, p, kappa, beta, mu, muz, sigma, sigmaz)
%calculate the exponent alpha
a = alpha(beta, mu, muz, sigma, sigmaz);
K = (mu/sigma)^2/2;
C = beta-K*p/(1-p);
q = 1/(p-1);
%equation for g(y, z) = w
F = @(w) ((1-p)/C)*(kappa^q - (w^2)^q)+y+z-z*(w^2/kappa)^(a-1);
w0 = sqrt(kappa/2);
[g,error] = fsolve(F, w0, optimset('Display','none'));
g = g^2;
end
function c = alpha(beta, mu, muz, sigma, sigmaz)
K = (mu^2/sigma^2)*0.5;
M = (sigmaz*mu)/sigma;
c = (M+K-beta)/(2*K)+sqrt((beta-M-K)^2+4*K*(beta-muz))/(2*K);
end