How to make my code run faster

Hi, i'm trying to solve a non linear system of five equations, they are very long equations that includes three parameter, lambda, sigma, and V, all three of the equation are of the form : 1/ lambda + 1/sigma + 1/V + 1/(lambda*sigma) + 1/(lambda^2) = C.

I tried solving it symbolically but it's too complicated, so i'm now trying to solve the sustem numerically. However when using the function solve, the code runs for almost 5 to 6 minutes then shut down. I have run my code step by step and it is the solve function that takes forever. I also tried using vpasolve and fsolve without succes, in those case the code doesn't work it keeps saying the object put as equatiosn are not symbolic.

Here's my code :

clc; close all; clear
%% Les données (input)
l = 30.6e-3;   % Epaisseur
f = 100e3;     % Frecq centrale
c0 = 340;      % Célérité ds l'air
rho_0 = 1.204; % densité air
eta = 18.5e-6; % Viscosité air
gamma = 1.4;   % constante GP
Pr = 0.707;    % Nombre de Prandlt
por = 0.95;    % Porosité
tor = 1.27;    % Tortuosité
%% Coefficients de calcul
xValue = double((2*pi*f*l*sqrt(2)/(2*c0))*(sqrt(eta/(rho_0*2*pi*f)))*(1+(gamma-1)/(3*sqrt(Pr))));
y1Value = double((l*eta/(rho_0*c0))*((gamma-1)/(3*Pr)+3));
y2Value = double((l*eta/(rho_0*c0))*(4*(gamma-1)/(3*sqrt(Pr))));
y3Value = double(-(l*eta/(2*rho_0*c0))*(1+(gamma-1)/(3*sqrt(Pr)))^2);
z1Value =  double((2*pi*f*l*sqrt(2)/(2*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(15*(gamma-1)/(4*27*Pr)+(15/4)));
z2Value =  double((2*pi*f*l*sqrt(2)/(2*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(6*(gamma-1)/(9*Pr)+6*(gamma-1)/(3*sqrt(Pr))));
z3Value = double(- (2*pi*f*l*sqrt(2)/(24*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))*((gamma-1)/(3*sqrt(Pr))+3));
z4Value = double(- (2*pi*f*l*sqrt(2)/(24*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))*(4*(gamma-1)/(3*sqrt(Pr))));
z5Value = double((2*pi*f*l*sqrt(2)/(8*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))^3);
z6Value = double(- (2*pi*f*l*sqrt(2)/(48*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))*((gamma-1)/(3*sqrt(Pr))+3));
z7Value = double(- (2*pi*f*l*sqrt(2)/(48*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))*(4*(gamma-1)/(3*sqrt(Pr))));
%% Calcul des paramètres
angles = [0, 30, 60];                   % Vecteur des angles d'incidence
coeff_T = [0.04501 0.02442 0.00252];    % Vecteur des coeff de transmission
% Initialisation des résultats
% parametres = zeros(3, length(angles));  % Une matrice pour stocker les 5 paramètres pour chaque angle
syms lambda sigma V;
parametres = [lambda, sigma, V];  % Une matrice pour stocker les 5 paramètres pour chaque angle
equations = sym(zeros(1, length(angles)));
% Boucle pour calculer les paramètres
for i = 1:length(angles)
    theta = angles(i);      % Angle d'incidence actuel
    T = coeff_T(i);         % Coefficient de transmission actuel
    E0 = por*sqrt(tor - sind(theta)^2)/(cosd(theta)*tor);
    sqrtTerm = sqrt(tor - sind(theta)^2);
    logTerm = log(abs(1 + E0)^2) + log(T) - log(4) - log(abs(E0));
    % Calcul des paramètres
    % Equation 1
    equations(i) = xValue*tor/(parametres(1)*sqrtTerm)...
        + y1Value*tor/(parametres(2)*sqrtTerm) + z1Value*tor/(parametres(3)*sqrtTerm)...
        + 1/(parametres(1))^2*( y2Value*tor/(sqrtTerm) + y3Value*tor^2/(sqrtTerm)^3)...
        + 1/(parametres(1)*parametres(2))*( z2Value*tor/(sqrtTerm) + z3Value*tor^2/(sqrtTerm)^3 + z6Value*tor^2/(sqrtTerm)^5)...
        +1/(parametres(1)^3)*( z5Value*tor^3/(sqrtTerm)^5 + z4Value*tor^2/(sqrtTerm)^3 + z7Value*tor^2/(sqrtTerm)^5) == logTerm ;
end
eqns_simpl = sym(zeros(1, length(angles)));
for i = 1:length(equations)
    eqns_simpl(i) = simplify(equations(i));
end
% Affichage des résultats
disp('Paramètres calculés :');
disp(parametres);
disp('Equations :');
disp(equations)
for i = 1:length(eqns_simpl)
    disp(class(eqns_simpl(i))); % Vérifiez si chaque élément est bien un 'sym'
end
assume(lambda > 0 & sigma > 0 & V > 0);
solution = vpasolve(equations, parametres, [1, 1, 1]);

Answer 1

[Torsten]

Maybe your system does not have a solution.

sol = fsolve(@fun,[1 1 0.1]);

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.

fun(sol)

ans = 1×3

-0.2988 -0.1292 0.2439

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

sol = 1./sol

sol = 1×3

-0.0001 2.1860 10.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function equations = fun(parametres)
%% Les données (input)
l = 30.6e-3;   % Epaisseur
f = 100e3;     % Frecq centrale
c0 = 340;      % Célérité ds l'air
rho_0 = 1.204; % densité air
eta = 18.5e-6; % Viscosité air
gamma = 1.4;   % constante GP
Pr = 0.707;    % Nombre de Prandlt
por = 0.95;    % Porosité
tor = 1.27;    % Tortuosité
%% Coefficients de calcul
xValue = (2*pi*f*l*sqrt(2)/(2*c0))*(sqrt(eta/(rho_0*2*pi*f)))*(1+(gamma-1)/(3*sqrt(Pr)));
y1Value = (l*eta/(rho_0*c0))*((gamma-1)/(3*Pr)+3);
y2Value = (l*eta/(rho_0*c0))*(4*(gamma-1)/(3*sqrt(Pr)));
y3Value = -(l*eta/(2*rho_0*c0))*(1+(gamma-1)/(3*sqrt(Pr)))^2;
z1Value =  (2*pi*f*l*sqrt(2)/(2*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(15*(gamma-1)/(4*27*Pr)+(15/4));
z2Value =  (2*pi*f*l*sqrt(2)/(2*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(6*(gamma-1)/(9*Pr)+6*(gamma-1)/(3*sqrt(Pr)));
z3Value = - (2*pi*f*l*sqrt(2)/(24*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))*((gamma-1)/(3*sqrt(Pr))+3);
z4Value = - (2*pi*f*l*sqrt(2)/(24*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))*(4*(gamma-1)/(3*sqrt(Pr)));
z5Value = (2*pi*f*l*sqrt(2)/(8*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))^3;
z6Value = - (2*pi*f*l*sqrt(2)/(48*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))*((gamma-1)/(3*sqrt(Pr))+3);
z7Value = - (2*pi*f*l*sqrt(2)/(48*c0))*((sqrt(eta/(rho_0*2*pi*f)))^3)*(1+(gamma-1)/(3*sqrt(Pr)))*(4*(gamma-1)/(3*sqrt(Pr)));
%% Calcul des paramètres
angles = [0, 30, 60];                   % Vecteur des angles d'incidence
coeff_T = [0.04501 0.02442 0.00252];    % Vecteur des coeff de transmission
% Initialisation des résultats
% parametres = zeros(3, length(angles));  % Une matrice pour stocker les 5 paramètres pour chaque angle
%syms lambda sigma V;
%parametres = [lambda, sigma, V];  % Une matrice pour stocker les 5 paramètres pour chaque angle
equations = zeros(1, length(angles));
% Boucle pour calculer les paramètres
for i = 1:length(angles)
    theta = angles(i);      % Angle d'incidence actuel
    T = coeff_T(i);         % Coefficient de transmission actuel
    E0 = por*sqrt(tor - sind(theta)^2)/(cosd(theta)*tor);
    sqrtTerm = sqrt(tor - sind(theta)^2);
    logTerm = log(abs(1 + E0)^2) + log(T) - log(4) - log(abs(E0));
    % Calcul des paramètres
    % Equation 1
    equations(i) = xValue*tor*parametres(1)/sqrtTerm...
        + y1Value*tor*parametres(2)/sqrtTerm + z1Value*tor*parametres(3)/sqrtTerm...
        + parametres(1)^2*( y2Value*tor/(sqrtTerm) + y3Value*tor^2/(sqrtTerm)^3)...
        + parametres(1)*parametres(2)*( z2Value*tor/(sqrtTerm) + z3Value*tor^2/(sqrtTerm)^3 + z6Value*tor^2/(sqrtTerm)^5)...
        +parametres(1)^3*( z5Value*tor^3/(sqrtTerm)^5 + z4Value*tor^2/(sqrtTerm)^3 + z7Value*tor^2/(sqrtTerm)^5) - logTerm ;
end
end

[Romain]

It should technically have solutions, but maybe there is a problem in the computing of my equation, because this is the computing of physical values.

Thank you for your help tho, i will try and go backto check if my equations are correct.

[Torsten]

I rewrote your equations by multiplying instead of dividing by your parameters. Remember to invert the obtained parameters as I did after receiving a solution from "fsolve".