Numerically Solving Complex Equation

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?

[Walter Roberson]

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]

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

[John]

@Walter Roberson thank you for a more optimized approach and @Torsten for a fully working, implementation of the equation. Since it seems that it is not possible to solve with Matlab, I will try it in Mathematica.

[Torsten]

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]

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]

1. 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.
2. Use vpasolve with Initial Guesses:Providing initial guesses close to the expected solutions can help vpasolve find a solution.
3. Adjust Tolerances if Necessary:If the solution is sensitive, you may need to adjust the tolerances for the solver.
4. 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.