Numerically Solving Complex Equation

2024 1 11
1 个回答
更新时间：2024 1 18
12 次查看（30 天）

0 个投票

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
Walter Roberson 2024-1-11
Ran in:
Ran in:
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
Torsten 2024-1-13
编辑：Torsten 2024-1-13
Ran in:
Ran in:
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
John 2024-1-16
@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
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
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.

请先登录，再进行评论。

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

回答（1 个）

Hassaan
Hassaan 2024-1-11

0 个投票

  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.

类别

在 帮助中心 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

产品

版本

R2022b

标签

提问:

John
John
2024-1-11

评论：

John
John
2024-1-18

Community Treasure Hunt

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

Start Hunting!

Translated by