Is this a correct way to use fsolve?

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

Answer 1

[Torsten]

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

[Salvatore Bianco]

Thanks, I solved by using lsqnonlin since I know the boundary where w should be.

[Salvatore Bianco]

Does this do what my previous script did? Since my function takes values only in (0, kappa], I use this so that i can specify does boundaries. It eliminated the problem of small residual complex parts, and it seems to work, but is this conceptually correct?

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

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("lsqnonlin", 'SpecifyObjectiveGradient', true);
fun = {F, J};
w0 = kappa/2;
g = lsqnonlin(fun, w0, 0, kappa, options);
end

[Salvatore Bianco]

I tried your method with w^2 but it gave me complex results again...

[Torsten]

I tried your method with w^2 but it gave me complex results again...

Can you include your code with an example for which this happens ? And don't use the "SpecifyObjectiveGradient" option (reason given below).

The point is that w^x will give complex values as soon as w becomes negative. In my opinion, this should be prevented by using w^2 instead of w in the function definition.

[Salvatore Bianco]

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)= g(y(j), z(i), p, kappa, beta, mu, muz, sigma, sigmaz);
end
end

function g = 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^q)+y+z-z*(w/kappa)^(a-1);
w0 = sqrt(kappa/2);
g = fsolve(F, w0);
g = sqrt(g);
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

[Torsten]

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

[Salvatore Bianco]

Sorry, I just edited it. I did not use g^2 since I am now solving for q^2. So when I do find it, I take the square root so I have g?

[Salvatore Bianco]

Sorry, I think the summer heat is hitting my brain, now I got it.