Clear Filters
Clear Filters

Fsolve Recommending Algorithm I Already Specified

3 views (last 30 days)
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.

Accepted Answer

Torsten
Torsten on 7 Dec 2023
Edited: Torsten on 7 Dec 2023
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 Comments
Matt J
Matt J on 7 Dec 2023
Edited: Matt J on 7 Dec 2023
So, do instead,
options = optimoptions(options, 'TolFun', 1e-6, 'TolX', 1e-6);
omega_gamma_solution = fsolve(____, options);

Sign in to comment.

More Answers (0)

Categories

Find more on Condensed Matter & Materials Physics in Help Center and File Exchange

Tags

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!