Gaps in contour plot using custom secant method (when solving for two-plume merger)

Hello all,

I'm working on a MATLAB script that computes the merging contour between two turbulent plumes, using a custom root-finding function called mnhf_secant. This function implements a supervisor-provided secant method to solve a nonlinear equation that describes the interface between the plumes. The governing equation being solved is located at the bottom of the code and is derived from plume interaction theory.

The issue arises in the else branch of the loop that iterates over polar angles (θ). When the code enters this branch, the contour plot generated shows visible gaps between the red line (plume boundary) and the black line (reference contour), indicating that the solution is not continuous or fails to resolve in those regions.

This breakdown does not occur in the first half of the angular domain (handled by the if condition), where the contours are computed correctly as seen by the continuous black line. The discontinuity seems localized to the else loop, despite using the same root-finding logic.I tried different guesses for the secant function but is not helping to solve for the gaps.

My main question is:

Is there a modification I could implement in the else branch to ensure a continuous contour across the full angular domain? I have attached the code below.

close all; clc; clf
global r0 k theta 
%% Parameter input, plume source of finite size
r0 = 1/3;
k_array = [0.3 0.5 0.7 1-r0^2 1]

k_array = 1×5

0.3000 0.5000 0.7000 0.8889 1.0000

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

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

Nt = 20001; % Must be odd for symmetry
theta_array = linspace(-pi/2, pi/2, Nt);
figure(1); hold on; box on
for jj = 1:length(k_array)
    k = k_array(jj);
    if k >= 1 - r0^2
        rho_array1 = zeros(1,length(theta_array));
        for ii = 1:length(theta_array)
            theta = theta_array(ii);
            rho_array1(ii) = mnhf_secant(@TwoPlumes_Equation, [1.5 2], 1e-6, 0);
            if rho_array1(ii)<0
                rho_array1(ii) = NaN;
            end
        end
        [x1, y1] = pol2cart(theta_array, rho_array1);
        plot(x1, y1, '-k'); axis square; xlim([0 2]); ylim([-1 1]);
    else
        rho_array1 = zeros(1,length(theta_array));
        for ii = 1:length(theta_array)
            theta = theta_array(ii);
            rho_array1(ii) = mnhf_secant(@TwoPlumes_Equation, [1.5 2], 1e-6, 0);
            if rho_array1(ii)<0
                rho_array1(ii) = NaN;
            end
        end
        [x1, y1] = pol2cart(theta_array, rho_array1);
        plot(x1, y1, '-k'); axis square; xlim([0 2]); ylim([-1 1]);
        
        % Identify breakdown angle where rho fails
        ij = find(isnan(rho_array1(1:round(end/2)))); 
        if ~isempty(ij)
            theta_min = theta_array(max(ij)) % closest valid angle before breakdown
            theta_array2 = linspace(theta_min, -theta_min, Nt); % symmetric gap domain
            rho_array2 = zeros(1, length(theta_array2));
            
            % Solve numerically in gap region
            for ii = 1:length(theta_array2)
                theta = theta_array2(ii);
                rho_array2(ii) = mnhf_secant(@TwoPlumes_Equation, [0.1 0.4], 1e-6, 0);
                if rho_array2(ii)<0
                    rho_array2(ii) = NaN;
                end
            end
            [x2, y2] = pol2cart(theta_array2, rho_array2);
            plot(x2, y2, 'r.'); axis square; xlim([0 2]); ylim([-1 1]);
        end
    end
    pos1 = [1 - r0, -r0, 2*r0, 2*r0];
    rectangle('Position', pos1, 'Curvature', [1 1], 'FaceColor', [0 0 0]);
end

theta_min = -0.3795

theta_min = -0.4486

theta_min = -0.5595

function r=mnhf_secant(Fun,x,tol,trace)
%MNHF_SECANT finds the root of "Fun" using secant scheme.
%
% Fun - name of the external function
% x - vector of length 2, (initial guesses)
% tol - tolerance
% trace - print intermediate results
%
% Usage  mnhf_secant(@poly1,[-0.5 2.0],1e-8,1)
%        poly1 is the name of the external function.
%        [-0.5 2.0] are the initial guesses for the root.
% Check inputs.
if nargin < 4, trace = 1;  end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 2)
    error('Please provide two initial guesses')
end
f = feval(Fun,x); % Fun is assumed to accept a vector
for i = 1:100
    x3 = x(1)-f(1)*(x(2)-x(1))/(f(2)-f(1));  % Update the guess.
    f3 = feval(Fun,x3);  % Function evaluation.
    if trace, fprintf(1,'%3i %12.5f %12.5f\n', i,x3,f3); end
    if abs(f3) < tol % Check for convergenece.
        r = x3; 
        return 
    else % Reset values for x(1), x(2), f(1) and f(2).
        x(1) = x(2); f(1) = f(2); x(2) = x3; f(2) = f3;
    end
end
r = NaN;
end
%% Two-Plume Interaction Equation (B3 from LF20B)
function f = TwoPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 - 2*rho.*cos(theta) + 1 - r0^2) .* ...
    (rho.^2 - 2*rho.*cos(theta + pi) + 1 - r0^2) - k^2;
end

[dpb]

Your problem is in what you don't provide, the solver code, mnhf_secant(). Nothing anybody here could possibly do without being able to run your code.

You need to set a breakpoint where it fails and look at the form of the equations in that area given the specific constants. It's likely there simply isn't a real, positive solution after some point.

r0 = 1/3;
k_array = [0.3 0.5 0.7 1-r0^2 1]

k_array = 1×5

0.3000 0.5000 0.7000 0.8889 1.0000

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

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

k_array >= (1 - r0^2)

ans = 1×5 logical array
   0   0   0   1   1

[Tashmid]

I attached the code mnhf_secant.

Answer 1

[Torsten]

The equation for rho is a polynomial of degree 4.

The MATLAB function "roots" solves for the four zeros of this equation.

I suspect that your function "mnhf_secant" in some cases either cannot find the correct root or there aren't any real-valued roots since all solutions are complex-valued.

Test it by calling "solve_TwoPlumes_equation" from below instead of "mnhf_secant". Note that r is an array with four elements (all the roots of the equation). So you will first have to check and sort out which of the four solutions (if any) makes sense physically and return this single specific value as "rho".

function rho = solve_TwoPlumes_Equation
  global r0 k theta
  p = [1, -2*(cos(theta+pi)+cos(theta)), 2*(1-r0^2)+4*cos(theta)*cos(theta+pi), -2*(1-r0^2)*(cos(theta)+cos(theta+pi)), (1-r0^2)^2-k^2];
  r = roots(p);
  rho = r(?);
end

[Walter Roberson]

After

r = roots(p);

you might want

r = r(imag(r)==0);

to gather only the real roots.

However, after that you run into the possibility that r will be empty (if there were no real roots) or that r might have 2 or 4 values. You might want

r = min(r, ComparisonMethod="abs");

for example, to get the r that has the smallest absolute value. Or you might want to follow with

tr = r(r>0);
if ~isempty(tr)
    r = tr(1);
else
    r = max(r);
end

to get the smallest positive root if one exists and otherwise the largest of the negative roots... or empty if there were no real roots.

[Tashmid]

Thanks @Walter Roberson for your response. I will try that and let you know.