How can I solve a algebraic equation with "fzero" ?

Solve "w" from an agebraic equation f [w; k,m,S,a]=0 . (All other symbols are constant values. The equation is shown in the end of the code)

(1) Question: How to use "fzero" to solve this equation numerically? (The code is shown below, but it doesn't work.)

(2) I've already know the solution of "w" , (shown in the figure below) "w" is a complex number, with real part "wr" and imaginary part "wi". I understand there are a lot solutions exist, but I would like to write a program and get the solution the same as the figure.

Thank you so much! I really appreciate your help.

Here is the program, which doesn't work:

clc
clear
%The equation is shown at the end. f[w; k,m,S,a]=0; where "k,m,S,a" are constant values. I would like to solve "w".
%define constant values in the equation
S=0.7; a=0; m=1;
k1=0.01:0.01:2.01;% k= 0 -> 2 with 0.01 interval.

for j=1:201
    k=k1(j); %set "k" value for each loop.
    w0=0;    %set initial point w0=0. The solution I want is positive complex values: w=wr + i*wi;
    eqn=@(w)f(w,k,m,S,a); 
    sol=fzero(eqn,w0);    %solve equation
    
end

% Define the function: solve "w" from f(w,k,m,S,a) = 0;
function y=f(w,k,m,S,a)
y=((w-m*S-(1+a)*k)+k)^2 ...
        *(-2*m*S+(w-m*S-(1+a)*k)*(k*sqrt((4*S^2-(w-m*S-(1+a)*k)^2)/(w-m*S-(1+a)*k)^2)) ...
        *1/2*(besselj(m-1,(k*sqrt((4*S^2-(w-m*S-(1+a)*k)^2)/(w-m*S-(1+a)*k)^2)))-besselj(m+1,(k*sqrt((4*S^2-(w-m*S-(1+a)*k)^2)/(w-m*S-(1+a)*k)^2))))...
        /besselj(m,(k*sqrt((4*S^2-(w-m*S-(1+a)*k)^2)/(w-m*S-(1+a)*k)^2))))...
        ...
        +(k*sqrt((4*S^2-(w-m*S-(1+a)*k)^2)/(w-m*S-(1+a)*k)^2))^2*(w-m*S-(1+a)*k)^3/k...
        *(-1/2)*(besselk(m-1,k)+besselk(m+1,k))/besselk(m,k);
end

[darova]

I tried this code

S=0.7; a=0; m=1;
k1=0.01:0.05:2.01;% k= 0 -> 2 with 0.01 interval.
w1 = k1;
[W,K] = meshgrid(w1,k1);
Z = W*0;

for i = 1:size(W,1)
    for j = 1:size(W,2)        
        Z(i,j) = f(W(i,j),K(i,j),m,S,a);
    end
end
clf
contour(W,K,real(Z),'k')
surface(W,K,real(Z),'edgecolor','none')
axis vis3d

Solution

[Lu Zhao]

Hi Darova,

Thank you so much for helping me! The idea to plot Z is great! I just tried this code out.

Just one more question: "w" is a complex number,(w=w_r+i*w_i), and do you have any ideas how to plot Z with complex number "w" ?

(And sorry for the confusion in the question, I will edit it right now to make it clear.)

Thanks again. I really appreciate your help! : D

Answer 1

[darova]

So you want to have 3 independent variables: w, and k

Here is one way: use isosurface

S=0.7; a=0; m=1;
k1=0.01:0.05:2.01;% k= 0 -> 2 with 0.01 interval.
w1 = k1;
[W1,W2,K] = meshgrid(w1,w1,k1);
Z = W1*0;
for i = 1:size(W1,1)
    for j = 1:size(W1,2)        
        for k = 1:size(W1,3)
            ww1 = W1(i,j,k);
            ww2 = W2(i,j,k);
            kk = K(i,j,k);
            Z(i,j,k) = f(ww1+1i*ww2,kk,m,S,a);
        end
    end
end
% clf
cla
patch( isosurface(W1,W2,K,imag(Z),0),...
    'facecolor','r',...
    'edgecolor','none');
patch( isosurface(W1,W2,K,real(Z),0),...
    'facecolor','b',...
    'edgecolor','none');
xlabel('W real')
ylabel('W imiginary')
zlabel('K variable')
legend('imiginary roots','real roots','location','north')
light
axis vis3d

[Lu Zhao]

Thank you so much Darova, this is a great visualization! I really appreciate your help.

I spend whole afternoon to learn about "patch" and "isosurface". Just wondered, is there any way to extract the values for "wr" and "wi" which makes f = 0, something like the figure showed in the question. (as "w1" and "w2" in you program)

[darova]

Sure, here is an example

[x,y,z,v] = flow;
p = patch(isosurface(x,y,z,v,-3));
set(p,'FaceColor','red','EdgeColor','none');
light
x = get(p,'xdata');
y = get(p,'ydata');
z = get(p,'zdata');
figure
line(x,y,z,'color','b')

[Lu Zhao]

Thank you ! I really appreciate it.

[Lu Zhao]

Hi Darova, hope you are having a great day so far.

I'm sorry for bothering you again. Just wondered, is there any way to get the complex conjugate pair of "w1 + i * w2" when abs(z)=0.

(I could see: w1 for real(z)=0 , and also w2 for imag(z)=0 from the code, but I have trouble with geting the conjugate pair that makes both real(z) and imag(z) equals zero (or 0.001) same time.)

Or is there any way to get solution for "w" looks like the same in figure below?

Thank you again. I really appreciate your help : D

[darova]

• Just wondered, is there any way to get the complex conjugate pair of "w1 + i * w2" when abs(z)=0.

Sure, use abs(). But in this case isosurface will not show anything because it's a line

One line you can find manually - intersection when imiginary is zero:

Another intersection can be found using surfaceIntersection

[Lu Zhao]

Thank you so much as always. I finally understand how to read this 3d figure, thanks a lot to your detailed picture. I couldn't thank you more. : D

Have a great day!

[darova]

my plesure