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]
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
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.5fn’, 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;
endHello 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]
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]);

% 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
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.5fn’, 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 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]
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]);

Cart

Cart