Plot showing gaps in numerical solution
I’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_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]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_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;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 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.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
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;
endI’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_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]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_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;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 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.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
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 I’m working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear
%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 – r0^2)^2; % Conservative lower bound
kU = 1.05*(1 – r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray)); % Store valid roots for each k in sweep
k_crit = NaN; % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
k = k_testArray(jj);
rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
if rho > 0 && isreal(rho) && isfinite(rho)
rho_valid(jj) = rho;
if ~contactDetected
k_crit = k;
contactDetected = true;
end
end
end
fprintf(‘Detected critical k value (first contact with θ = π/4): k_crit = %.8fnn’, k_crit);
%% Parameter input: k values
Nt = 10001; % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]
%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5]; % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
k = k_array(jj);
% For merged Case
if k >= k_crit
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_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]);
% Plot θ = π/4 lines for visual reference
x_diag = linspace(0, 2, Nt); y_diag = x_diag;
plot(x_diag, y_diag, ‘r–‘); plot(x_diag, -y_diag, ‘r–‘);
% Area estimation using trapezoidal rule
iq = find(y1 >= y_diag); % intersection
if ~isempty(iq)
area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
area(jj) = area(jj) + x1(min(iq)-1)^2;
end
else
% Before Merging Case
rho_array1 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array(ii);
rho_array1(ii) = mnhf_secant(@FourPlumes_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;
% Detect breakdown region and recompute in "gap"
ij = find(isnan(rho_array1(1:round(end/2))));
if ~isempty(ij)
theta_min = theta_array(max(ij));
theta_array2 = linspace(theta_min, -theta_min, Nt);
rho_array2 = zeros(1, Nt);
for ii = 1:Nt
theta = theta_array2(ii);
rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 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.’);
ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
plot(x, y, ‘b–‘); plot(x, -y, ‘b–‘);
% Area from symmetric region
area(jj) = 2 * trapz(x, y);
end
end
% Plot plume source circles (n = 4)
pos1 = [1 – r0, -r0, 2*r0, 2*r0]; % (1, 0)
pos2 = [-1 – r0, -r0, 2*r0, 2*r0]; % (-1, 0)
pos3 = [-r0, 1 – r0, 2*r0, 2*r0]; % (0, 1)
pos4 = [-r0, -1 – r0, 2*r0, 2*r0]; % (0, -1)
rectangle(‘Position’, pos1, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos2, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos3, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
rectangle(‘Position’, pos4, ‘Curvature’, [1 1], ‘FaceColor’, ‘k’);
end
%% Function: Four-Plume Equation (Li & Flynn B3) —
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 – 2*rho*cos(theta) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi/2) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + pi) + 1 – r0^2) .* …
(rho.^2 – 2*rho*cos(theta + 1.5*pi) + 1 – r0^2) – k^2;
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 numerical solution, nonlinear-equations, root-finding, polar-coordinates, contour-gaps MATLAB Answers — New Questions
