Fsolve Recommending Algorithm I Already Specified

2 次查看（过去 30 天）
John
John 2023-12-7
评论： John 2023-12-8
I'm trying to solve the dispersion relation in plasma and with this code, for some reason, even when I specify that solver as "levenberg-marquardt", it still gives me a warning saying this:
> In fsolve (line 345)
In GrowthRateSolving (line 41)
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead.
Not sure why it's doing that? Code:
% Constants (adjust as needed)
k = 1e7; % wave number (1/m)
lambda_De = 1e-3; % Debye length for electrons (m)
lambda_Di = 1e-2; % Debye length for ions (m)
me = 9.10938356e-31; % mass of electron (kg)
mi = 1.6726219e-27; % mass of ion (proton, kg)
kappa = 1.380649e-23; % Boltzmann constant (J/K)
Ti = 1e4; % ion temperature (K), kept constant
Te_Ti_range = linspace(0.01, 25, 100); % Example definition, adjust as needed
Te_Ti_sample = 1.05877268798617;
gamma_over_omega_sample = 0.805171624713958;
% Estimate Te from the sample ratio
Te_sample = Te_Ti_sample * Ti;
% Initial guess for omega
omega_initial = 1e8; % Adjust based on your knowledge of the system
% Solve for omega
%options = optimoptions('fsolve', 'Display', 'none', 'TolFun', 1e-6, 'TolX', 1e-6);
%options = optimoptions('fsolve', 'Algorithm','levenberg-marquardt', 'Display', 'none', 'TolFun', 1e-6, 'TolX', 1e-6);
options = optimoptions('fsolve', 'Algorithm', 'levenberg-marquardt', 'Display', 'none');
%options = optimoptions('fsolve', 'Display', 'none');
% Solve for omega
omega_solution = fsolve(@(omega) omega_relation(omega, gamma_over_omega_sample, Te_sample, Ti, me, mi, kappa, k, lambda_De, lambda_Di), omega_initial, options);
% Use omega to find the initial guess for gamma
gamma_initial_guess = -gamma_over_omega_sample * omega_solution;
% Loop over Te/Ti
for i = 1:length(Te_Ti_range)
Te = Te_Ti_range(i) * Ti; % Adjust Te based on Te/Ti ratio
% Solve for omega and gamma using fsolve
omega_gamma_initial = [1e6, 1e3]; % Example initial guesses
options = optimoptions('fsolve', 'Display', 'none', 'TolFun', 1e-6, 'TolX', 1e-6);
omega_gamma_solution = fsolve(@(omega_gamma) dispersion_relation(omega_gamma, Te, Ti, me, mi, kappa, k, lambda_De, lambda_Di), omega_gamma_initial, options);
omega = omega_gamma_solution(1);
gamma = omega_gamma_solution(2);
% Calculate -gamma/omega and store it
minus_gamma_over_omega_results(i) = -gamma / omega;
end
% Plotting
loglog(Te_Ti_range, minus_gamma_over_omega_results);
xlabel('Te/Ti');
ylabel('-gamma/omega');
axis([0.01 25 0.01 1]); % Setting the axis limits
grid on;
Supporting FUnctions:
function Z = plasma_dispersion(zeta)
% Parameters for the contour integration
R = 1.5; % Radius of semicircle
N = 1000; % Number of points for numerical integration
% Constructing a contour that avoids the singularity at zeta
theta = linspace(0, pi, N);
z = R * exp(1i * theta);
z = z + real(zeta);
% Integration along the contour
integrand = @(z) exp(-z.^2) ./ (z - zeta);
Z = trapz(z, integrand(z)) / sqrt(pi);
end
function F = dispersion_relation(omega_gamma, Te, Ti, me, mi, kappa, k, lambda_De, lambda_Di)
omega = omega_gamma(1);
gamma = omega_gamma(2);
x_e = sqrt(me/(2*kappa*Te)) * omega / k;
y_e = sqrt(me/(2*kappa*Te)) * gamma / k;
zeta_e = x_e + 1i*y_e;
x_i = sqrt(mi/(2*kappa*Ti)) * omega / k;
y_i = sqrt(mi/(2*kappa*Ti)) * gamma / k;
zeta_i = x_i + 1i*y_i;
sum_s = (1/(k*lambda_De)^2) * 0.5 * real(plasma_dispersion(zeta_e)) + ...
(1/(k*lambda_Di)^2) * 0.5 * real(plasma_dispersion(zeta_i));
F = 1 - sum_s;
end
function F = omega_relation(omega, gamma_over_omega_sample, Te, Ti, me, mi, kappa, k, lambda_De, lambda_Di)
gamma = -gamma_over_omega_sample * omega;
% Calculate zeta for electrons
x_e = sqrt(me/(2*kappa*Te)) * omega / k;
y_e = sqrt(me/(2*kappa*Te)) * gamma / k;
zeta_e = x_e + 1i*y_e;
% Calculate zeta for ions
x_i = sqrt(mi/(2*kappa*Ti)) * omega / k;
y_i = sqrt(mi/(2*kappa*Ti)) * gamma / k;
zeta_i = x_i + 1i*y_i;
% Summation over species (electrons and ions)
sum_s = (1/(k*lambda_De)^2) * 0.5 * real(plasma_dispersion(zeta_e)) + ...
(1/(k*lambda_Di)^2) * 0.5 * real(plasma_dispersion(zeta_i));
% Dispersion relation equation
F = 1 - sum_s;
end
I've already tested the solver with a simply quadratic to make sure it is working and it does for that:
% Initial guess for the solution
x_initial = 2; % Adjust this to test different initial guesses
% Options for fsolve
options = optimoptions('fsolve', 'Display', 'iter', 'Algorithm', 'levenberg-marquardt');
% Solve the equation
x_solution = fsolve(@simplified_test_equation, x_initial, options);
% Display the solution
disp(['The solution is x = ', num2str(x_solution)]);
% Define the simplified test equation
function F = simplified_test_equation(x)
F = x^2 - 4*x + 3;
end
And my more complex code is running and functioning, so not sure what else could be the issue.

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

采纳的回答

Torsten
Torsten 2023-12-7
编辑：Torsten 2023-12-7
You overwrite the previously set options with this line
options = optimoptions('fsolve', 'Display', 'none', 'TolFun', 1e-6, 'TolX', 1e-6);
before calling "fsolve" in the loop.
  2 个评论
Matt J
Matt J 2023-12-7
编辑：Matt J 2023-12-7
So, do instead,
options = optimoptions(options, 'TolFun', 1e-6, 'TolX', 1e-6);
omega_gamma_solution = fsolve(____, options);

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by