Numerically Solving Complex Equation

Question

John 2024-1-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2068986-numerically-solving-complex-equation

评论： John 2024-1-18

I am trying to solve the following equation numerically:

Where I am trying to solve for omega in terms of k and both are complex and nonzero. Lamda_D and v_Ts are both real. Both contain real and imaginary components. Z is the plasma dispersion function. Here is how I implemented it in code:

% User Inputs
ne = (1E+8); %charge #/m^3
ni = (1E+8); %charge #/m^3
percIon = 1; % Percentage of molecules that are ions
d = 70; % Electrode spacing in (cm)
z = 1; % Effective charge of ions in plasma
Ti = 465;
Te = 471;
% Natural Constants
u0 = pi*4E-7; %(H/m)
ep0 = 8.8542E-12; %Vacuum constant (Farads/meter)
q = 1.60218E-19; %Absolute charge of electron (Coulombs)
k_B = 1.380649E-23; %Joule/K
m_i = 1.667E-27;
m_e = 9.11E-31;
% Calculations 
v_Ti = sqrt((3*k_B*Ti)/m_i);
v_Te = sqrt((3*k_B*Te)/m_e);
%In (K)
Num = (ep0*k_B)/q^2;
Den = ne/Te + z^2*ni/Ti;
lambdaD = sqrt(Num/(percIon*Den));
syms omega k
squi_i = omega / (k * v_Ti * sqrt(2));
squi_e = omega / (k * v_Te * sqrt(2));
plasDispFun_Di = -2*(1 + squi_i*1i*sqrt(pi)*exp(-squi_i^2)*(1 + erf(1i*squi_i)));
plasDispFun_De = -2*(1 + squi_e*1i*sqrt(pi)*exp(-squi_e^2)*(1 + erf(1i*squi_e)));
eps = 1 - plasDispFun_De/(2*k^2*lambdaD) - plasDispFun_Di/(2*k^2*lambdaD);
sol = vpasolve(eps,[omega k]);

I tried using vpasolve but I get a solution that is empty.

[sol.omega sol.k]
ans =
Empty sym: 0-by-2

Does anyone know what I may be doing incorrectly here? Or is this not the proper function to use for this application?

5 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2024-1-11

在 MATLAB Online 中打开

A more precise calculation

Q = @(v) sym(v);

% User Inputs

ne = Q(1E+8); %charge #/m^3

ni = Q(1E+8); %charge #/m^3

percIon = Q(1); % Percentage of molecules that are ions

d = Q(70); % Electrode spacing in (cm)

z = Q(1); % Effective charge of ions in plasma

Ti = Q(465);

Te = Q(471);

% Natural Constants

u0 = Q(pi)*Q(4E-7); %(H/m)

ep0 = Q(88542)*Q(10)^-16; %Vacuum constant (Farads/meter)

q = Q(160218)*Q(10)^-24; %Absolute charge of electron (Coulombs)

k_B = Q(1380649)*Q(10)^-29; %Joule/K

m_i = Q(1667)*Q(10)^-30;

m_e = Q(911)*Q(10)^-33;

% Calculations

v_Ti = sqrt((Q(3)*k_B*Ti)/m_i);

v_Te = sqrt((Q(3)*k_B*Te)/m_e);

%In (K)

Num = (ep0*k_B)/q^2;

Den = ne/Te + z^2*ni/Ti;

lambdaD = sqrt(Num/(percIon*Den));

syms omega k

squi_i = omega / (k * v_Ti * sqrt(2));

squi_e = omega / (k * v_Te * sqrt(2));

plasDispFun_Di = -2*(1 + squi_i*1i*sqrt(pi)*exp(-squi_i^2)*(1 + erf(1i*squi_i)));

plasDispFun_De = -2*(1 + squi_e*1i*sqrt(pi)*exp(-squi_e^2)*(1 + erf(1i*squi_e)));

eps = simplify(1 - plasDispFun_De/(2*k^2*lambdaD) - plasDispFun_Di/(2*k^2*lambdaD))

eps =

symvar(eps)

ans =

sol = vpasolve([real(eps), imag(eps)] ,[omega k])

sol = struct with fields:

omega: [0×1 sym] k: [0×1 sym]

Unfortunately, the numeric erf function cannot handle complex values

fun = matlabFunction([real(eps), imag(eps)], 'Vars', {[omega k]});
soln = fsolve(fun, [.1, .2])
Error using erf
Input must be real and full.

Error in symengine>@(in1)[real(2.931550633783428e+7.*1.0./in1(:,2).^2.*((2.671100866178587e+5.*in1(:,1).*exp(1.0./in1(:,2).^2.*in1(:,1).^2.*(-4.327609982164354e-8)).*(erf((2.671100866178587e+5.*in1(:,1).*7.788140339827871e-10i)./in1(:,2))+1.0).*2.760823867344101e-9i)./in1(:,2)+2.0)).*1.615729297377441e-7+real(2.931550633783428e+7.*1.0./in1(:,2).^2.*((6.284426537449539e+5.*in1(:,1).*exp(1.0./in1(:,2).^2.*in1(:,1).^2.*(-2.334871245616549e-11)).*(erf((6.284426537449539e+5.*in1(:,1).*7.688928360976561e-12i)./in1(:,2))+1.0).*2.725654136549908e-11i)./in1(:,2)+2.0)).*1.615729297377441e-7+1.0,imag(2.931550633783428e+7.*1.0./in1(:,2).^2.*((2.671100866178587e+5.*in1(:,1).*exp(1.0./in1(:,2).^2.*in1(:,1).^2.*(-4.327609982164354e-8)).*(erf((2.671100866178587e+5.*in1(:,1).*7.788140339827871e-10i)./in1(:,2))+1.0).*2.760823867344101e-9i)./in1(:,2)+2.0)).*1.615729297377441e-7+imag(2.931550633783428e+7.*1.0./in1(:,2).^2.*((6.284426537449539e+5.*in1(:,1).*exp(1.0./in1(:,2).^2.*in1(:,1).^2.*(-2.334871245616549e-11)).*(erf((6.284426537449539e+5.*in1(:,1).*7.688928360976561e-12i)./in1(:,2))+1.0).*2.725654136549908e-11i)./in1(:,2)+2.0)).*1.615729297377441e-7]

Error in fsolve (line 270)
            fuser = feval(funfcn{3},x,varargin{:});

Caused by:
    Failure in initial objective function evaluation. FSOLVE cannot continue.

Torsten 2024-1-13

编辑：Torsten 2024-1-13

在 MATLAB Online 中打开

Use erf(i*x) = i*erfi(x) for real x. But MATLAB cannot find a solution:

% User Inputs
ne = (1E+8); %charge #/m^3
ni = (1E+8); %charge #/m^3
percIon = 1; % Percentage of molecules that are ions
d = 70; % Electrode spacing in (cm)
z = 1; % Effective charge of ions in plasma
Ti = 465;
Te = 471;
% Natural Constants
u0 = pi*4E-7; %(H/m)
ep0 = 8.8542E-12; %Vacuum constant (Farads/meter)
q = 1.60218E-19; %Absolute charge of electron (Coulombs)
k_B = 1.380649E-23; %Joule/K
m_i = 1.667E-27;
m_e = 9.11E-31;
% Calculations 
v_Ti = sqrt((3*k_B*Ti)/m_i);
v_Te = sqrt((3*k_B*Te)/m_e);
%In (K)
Num = (ep0*k_B)/q^2;
Den = ne/Te + z^2*ni/Ti;
lambdaD = sqrt(Num/(percIon*Den));
squi_i = @(k,omega)omega / (k * v_Ti * sqrt(2));
squi_e = @(k,omega)omega / (k * v_Te * sqrt(2));
plasDispFun_Di = @(k,omega)-2*(1 + squi_i(k,omega)*1i*sqrt(pi)*exp(-squi_i(k,omega)^2)*(1 + 1i*erfi(squi_i(k,omega))));
plasDispFun_De = @(k,omega)-2*(1 + squi_e(k,omega)*1i*sqrt(pi)*exp(-squi_e(k,omega)^2)*(1 + 1i*erfi(squi_e(k,omega))));
eps = @(k,omega)1 - plasDispFun_De(k,omega)/(2*k^2*lambdaD) - plasDispFun_Di(k,omega)/(2*k^2*lambdaD);
sol = fsolve(@(x)[real(eps(x(1),x(2))),imag(eps(x(1),x(2)))],[.1 .2])
No solution found.

fsolve stopped because the problem appears regular as measured by the gradient,
but the vector of function values is not near zero as measured by the
value of the function tolerance.
sol = 1×2
  351.3438    3.1891

Torsten 2024-1-16

If you make a surface plot of abs(eps) over an arbitrary 2d-domain, you will see that abs(eps) does not have a zero. So it's not due to MATLAB, but due to your function definition that "fsolve" does not succeed.

John 2024-1-18

That's really odd. I also plotted it to double check. This is probably the simplest form of the dielectric function which has been solved many times before. I may have written the equation incorrectly and have missed something. Thanks for pointing that out.

请先登录，再进行评论。

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

Answer 1

Hassaan 2024-1-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2068986-numerically-solving-complex-equation#answer_1388191

在 MATLAB Online 中打开

Define the Equations Symbolically:Use symbolic math to define your equations. Since you have complex numbers, make sure to use syms with the 'complex' option.
Use vpasolve with Initial Guesses:Providing initial guesses close to the expected solutions can help vpasolve find a solution.
Adjust Tolerances if Necessary:If the solution is sensitive, you may need to adjust the tolerances for the solver.
Separate Real and Imaginary Parts:If vpasolve has trouble with the complex equation, you might need to separate the real and imaginary parts of the equation.

% Define the symbolic variables and equations
syms k real;
syms omega complex;
% Given constants (you will need to define these with actual values)
lambda_D = ...; % replace with actual value
v_Ts = ...; % replace with actual value
% Define xi_s in terms of omega and k
xi_s = omega / (k * v_Ts * sqrt(2));
% Define Z'(xi_s)
Z_prime = -2 * (1 + xi_s * i * sqrt(pi) * exp(-xi_s^2) * (1 + erf(i * xi_s)));
% Main equation to solve
eqn = 1 - (Z_prime / (2 * k^2 * lambda_D^2)) == 0;
% Solve the equation numerically for omega
% Provide an initial guess for omega if known
init_guess = ...; % replace with your initial guess for omega
options = optimset('Display', 'off', 'TolFun', 1e-6, 'TolX', 1e-6); % Adjust tolerances as needed
[omega_sol, ~, exitflag] = vpasolve(eqn, omega, init_guess, 'Options', options);
% Check if a solution was found
if exitflag
    fprintf('Solution found: omega = %s\n', vpa(omega_sol, 6));
else
    fprintf('No solution found.\n');
end

Replace the ... with the actual values you have for the constants and the initial guess for omega.

Please make sure to use your actual values for lambda_D, v_Ts, and an initial guess for omega that is close to where you expect the root to be. The initial guess can significantly affect the success of vpasolve in finding a solution. If you're unsure of the initial guess, you might need to explore the parameter space or use a numerical solver like fsolve that can handle complex numbers with appropriate options set.

---------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

Professional Interests

Technical Services and Consulting
Embedded Systems | Firmware Developement | Simulations
Electrical and Electronics Engineering

Feel free to contact me.

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

请先登录，再进行评论。

Numerically Solving Complex Equation

5 个评论
显示 3更早的评论隐藏 3更早的评论

回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Numerically Solving Complex Equation

5 个评论 显示 3更早的评论隐藏 3更早的评论

回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论

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