How do I make my code faster (currently sym.vpasolve takes 43,434 seconds for 30,000 calls)

% Universal Constants
F = 96485; % farads
R = 8.31; % gas
NA = 6.02e23; % avagadros
T = 298; % in kelvin
% Cell constants
Length = 2.95 * 1e-6;
Width = 1.07 * 1e-6;
Radius = Width / 2;
% Definition of surface and volume (for a spherocylinder) = π•r²•(l-2r) [volume of cylindrical section whose length is l-2r] + 4/3•π•r³ [volume of speherical section]
S = pi * Width * Length;
Vol = pi * (Width / 2)^2 * (Length - Width / 3);
% Membrane specific capacitance
Cm = 6.5e-3;
constant = (F * Vol) / (S * Cm); % FV/SCm constant in memb potential eqn
eta = F / (R * T); % F/RT
pKa_w = 14; % pKa of water
% External environments
pHe = 7; % Neutral environment
% Concentration of protons [mM]
H_c_e = 10^(pHe) * 1e3;
% Concentration of hydroxides [mM]
OH_c_e = 10^(pHe - pKa_w) * 1e3;
% intracellular trapped ions conc in mM
z_y = [-10, 0, 10];
% In the case of pHe= 7= pHi, PMF= psi (membrane voltage)
pmf = linspace(-0.24, 0, 100); % setting 101 regularly spaced values for PMF
CA_e = linspace(0.1, 100, 100); %0.1 mM to 100mM in total contribution from [Na] and [Cl] ions 
total_iterations = length(CA_e) * length(z_y) * length(pmf);
count = 0;
% Preallocate cell arrays to store data
%variables_values = cell(total_iterations, 1);
%solutions = cell(total_iterations, 1);
syms x;
for CA_e_i = CA_e
    % 
    C_c_e = 0.5*CA_e_i; % in mM
    % fractional ionic composition
    if CA_e_i == 0
        continue;
    else
        alpha_c = C_c_e / CA_e_i;
    end
    % Concentration of anions
    A_c_e = 0.5*CA_e_i; % in mM
    if CA_e_i == 0
        continue;
    else
        alpha_a = A_c_e / CA_e_i;
    end
    for z_y_i = z_y
        if CA_e_i == 0
            continue;
        else
            alpha_y = z_y_i / CA_e_i;
        end
        for pmf_i = pmf
            psi = pmf_i;
            count = count + 1;
            % Define the equation, Eq is an equation object, first term is the rhs and the second term is the LHS (normally what it is equal to)
            if CA_e_i == 0
                continue;
            else
                lhs = (psi) / constant * CA_e_i; % from eqn S42 in SI
            end
            % assuming only cation exporting antiporters work
            equation = alpha_y + alpha_c * exp(eta * (x - psi) - alpha_a * exp(eta * psi)) == lhs;
            solution = vpasolve(equation);
            solution_double = double(solution);
            % Store variable values and solution
            variables_values{count} = struct('CA_e_i', CA_e_i, 'z_y_i', z_y_i, 'pmf_i', pmf_i);
            
            %this has empty values for places where the solution is not a
            %real number (since the equation represents a physical system)
            solutions{count} = solution_double;
        end
    end
end

Answer 1

[Ronit]

Hi Keshav,

I see that you are looking for ways to optimize MATLAB code that currently takes an excessively long time to execute, specifically due to repeated calls to ‘vpasolve’ within nested loops.

I can suggest some strategies to make your code run faster:

1. Reduce Precision When Possible: ‘vpasolve’ allows you to specify the precision of the solution. If you do not need very high precision, reducing it can significantly speed up the computation. Use the 'Digits' option to lower the precision. You can refer to the documentation: https://www.mathworks.com/help/symbolic/sym.vpasolve.html

solution = vpasolve(equation, x, 'Random', true, 'Digits', 4); % Example with reduced precision  

2. Preallocate memory for your variables outside the loop. It seems you have commented out the preallocation lines; ensure they are active and correctly sized to avoid dynamic resizing during the loop.

3. Profiling: Use MATLAB’s profiler (profile on; profile viewer) to identify the exact bottlenecks in your code.

4. Parallel Computing: You can use the Parallel Computing Toolbox to distribute the iterations across multiple cores. This is particularly effective for loop structures that are independent of each iteration. More about the Parallel Computing Toolbox can be gathered here: https://www.mathworks.com/products/parallel-computing.html

I hope these strategies will aid in making your code run faster.