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

1 次查看（过去 30 天）
Lu Zhao
Lu Zhao 2020-5-1
评论： darova 2020-5-4
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
  2 个评论
darova
darova 2020-5-1
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
Lu Zhao 2020-5-1
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

请先登录，再进行评论。

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

采纳的回答

darova
darova 2020-5-1
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
  7 个评论
显示 5更早的评论
Lu Zhao
Lu Zhao 2020-5-4
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!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by