Fsolve Recommending Algorithm I Already Specified

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.

Answer 1

[Torsten]

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.

[Matt J]

So, do instead,

    options = optimoptions(options, 'TolFun', 1e-6, 'TolX', 1e-6);
    omega_gamma_solution = fsolve(____, options);

[John]

Ah, haha thank you