Matching solutions of the heat equation in cylindrical coordinates

Difficulty Matching PDE Solutions at Using MATLAB Integration

Hello,

I am working on solving a boundary value problem involving a PDE (heat equation in cylindrical coordinates with two Neumann conditions and Dirac delta function as an IC) with different regimes for and .

where r is the radial coordinate, ϵ is of the order of 0.4 is of the order of 1.4 , is of the order of 2.4 :

κ is of the order of 1 and D in the order of 0.33. The BCs are: . The IC is a Dirac delta function at : .

The problem requires matching solutions at , using continuity condition for the function and demanding the integration equals to 1 of the eigenfunctions (normalization condition). However, my MATLAB implementation does not yield solutions that match well (there is a certain jump there) at .

I suspect there may be an issue in my approach and would appreciate any guidance. These are the solutions for each region :

The general solution inside :

where λ are the eigenvalues.

We have the BC:

.

where and are coefficients to be determined, and is the confluent hypergeometric function defined as

Continue:

Outside of :

where and are the zeroth order first and second kind Bessel function and A and B are coefficients to be determined and α is the square root of the eigenvalue λ as .

Boundary condition:

Substituting:

Final solution:

For matching we demand:

and the normalization condition as (when assuming $C_{2}$ to be zero which should be correct (as we care only about the ratio of the coefficients as later on they are ''swallowed'' when implementing the IC)):

So I start with isolating A:

Then, to obtain the eignevalues I substitute in the continuity condition:

Then I substitute what I got for A (which I'm not sure if the plus sign, minus sign or both is/are correct) and solve for the eigenvalues with the build-in MATLAB function . I can clrearly see from my plots that there is a relatively high jump between the two regions which I'm not happy about. I'll be glad if someone could advice from experience in matching solutions of PDEs.

Below is the MATLAB code. Thanks! Roi

close all
clear all
clc
% Parameters
R0 = 1.6;       % Inner radius
R_star = 2.4;   % Outer radius
D = 0.33;       % Diffusion coefficient
epsilon = 0.4;  % Small cutoff
kappa = 1;      % Reaction coefficient
% Special functions (placeholders for the confluent hypergeometric function M)
M = @(a, b, x) hypergeom(a, b, x); % MATLAB has the hypergeom function
J0 = @(x) besselj(0, x);
J1 = @(x) besselj(1, x);
Y0 = @(x) bessely(0, x);
Y1 = @(x) bessely(1, x);
% Define R1(r)^2 as a function of alpha
R1_squared = @(r, alpha) ( ...
    ( (1 - (alpha.^2 * D) / (2 * kappa)) / 2 * M(2 - (alpha.^2 * D) / (2 * kappa), 3, (kappa / (2 * D)) * (2 * epsilon)^2) ...
    ./ ((alpha.^2 * D) / (2 * kappa) * M(1 - (alpha.^2 * D) / (2 * kappa), 2, (kappa / (2 * D)) * (2 * epsilon)^2)) ...
    .* M(-(alpha.^2 * D) / (2 * kappa), 1, (kappa / (2 * D)) * r.^2) ...
    + M(1 - (alpha.^2 * D) / (2 * kappa), 2, (kappa / (2 * D)) * r.^2) ) ).^2;
% Define R2(r)^2 as a function of alpha
R2_squared = @(r, alpha) ( ...
    besselj(0, alpha .* r) - (besselj(1, alpha .* R_star) ./ bessely(1, alpha .* R_star)) .* bessely(0, alpha .* r) ).^2;
% Define the equation to solve for alpha
equation = @(alpha) ...
    ( ...
    besselj(0, alpha .* R0) - ...
    (besselj(1, alpha .* R_star) ./ bessely(1, alpha .* R_star)) .* bessely(0, alpha .* R0) ...
    ) .* sqrt( ...
    (1 - integral(@(r) R1_squared(r, alpha) .* exp(-kappa / (2 * D) * r.^2), 2 * epsilon, R0)) ./ ...
    integral(@(r) R2_squared(r, alpha) .* r, R0, R_star) ...
    ) ...
    - ...
    ( (1 - (alpha.^2 .* D) ./ (2 .* kappa)) ./ 2 ) .* ...
    M(2 - (alpha.^2 .* D) ./ (2 .* kappa), 3, (kappa ./ (2 .* D)) .* (2 .* epsilon).^2) ./ ...
    ((alpha.^2 .* D) ./ (2 .* kappa) .* M(1 - (alpha.^2 .* D) ./ (2 .* kappa), 2, (kappa ./ (2 .* D)) .* (2 .* epsilon).^2)) .* ...
    M(-(alpha.^2 .* D) ./ (2 .* kappa), 1, (kappa ./ (2 .* D)) .* R0.^2) + ...
    M(1 - (alpha.^2 .* D) ./ (2 .* kappa), 2, (kappa ./ (2 .* D)) .* R0.^2);

% Initial guesses for finding roots
% initial_guesses = [linspace(1e-4, 1e-3, 10),linspace(1e-3, 1e-2, 10),linspace(1e-2, 1e-1, 10),linspace(1e-1, 1e0, 10),linspace(1, 50, 100)]; % Array of initial guesses
initial_guesses = [linspace(1e-1, 1e0, 10),linspace(1, 50, 50)]; % Array of initial guesses
alpha_values = []; % Store distinct roots
tolerance = 1e-2;  % Minimum difference between roots
% Find eigenvalues (alpha values)
for guess = initial_guesses
    root = fsolve(equation, guess, optimoptions('fsolve', 'Display', 'none'));
    if isempty(alpha_values) || all(abs(alpha_values - root) > tolerance)
        alpha_values = [alpha_values, root];
    end
    if length(alpha_values) >= 10, break; end
end
% Double-check alpha values
valid_alpha_values = [];
for i = 1:length(alpha_values)
    if abs(equation(alpha_values(i))) < 1e-2
        valid_alpha_values = [valid_alpha_values, alpha_values(i)];
    end
end
disp('Validated alpha_i values (|equation(alpha)| < 1e-3):');
disp(valid_alpha_values);
% Radial functions
r_range1 = linspace(2 * epsilon, R0, 100); % Region 1
r_range2 = linspace(R0, R_star, 100);      % Region 2
% Define radial functions
R1 = @(r, alpha) ...
    ( (1 - (alpha^2 * D) / (2 * kappa)) / 2 ) .* ...
    M(2 - (alpha^2 * D) / (2 * kappa), 3, (kappa / (2 * D)) * (2 * epsilon)^2) ./ ...
    ((alpha^2 * D) / (2 * kappa) .* M(1 - (alpha^2 * D) / (2 * kappa), 2, (kappa / (2 * D)) * (2 * epsilon)^2)) .* ...
    M(-(alpha^2 * D) / (2 * kappa), 1, (kappa / (2 * D)) * r.^2) + ...
    M(1 - (alpha^2 * D) / (2 * kappa), 2, (kappa / (2 * D)) * r.^2);
R2 = @(r, alpha) ...
    besselj(0, alpha * r) - (besselj(1, alpha * R_star) / bessely(1, alpha * R_star)) * bessely(0, alpha * r);
% Compute integrals and normalize radial functions
coefficients1 = zeros(size(valid_alpha_values));
coefficients2 = zeros(size(valid_alpha_values));
% Progress bar for coefficients computation
h1 = waitbar(0, 'Calculating coefficients...');
for i = 1:length(valid_alpha_values)
    alpha = valid_alpha_values(i);
    % Compute R1 and R2 at R0
    R1_at_R0 = R1(R0, alpha);
    R2_at_R0 = R2(R0, alpha);
    % Compute the coefficients directly using the integral formulas
    coefficients1(i) = (1 / (2 * pi * R0)) * (R1_at_R0 * exp(-kappa ./ (2 .* D) .* R0.^2) / trapz(r_range1, R1(r_range1, alpha).^2 .* exp(-kappa ./ (2 .* D) .* r_range1.^2)));
    coefficients2(i) = (1 / (2 * pi)) * (R2_at_R0 / trapz(r_range2, R2(r_range2, alpha).^2 .* r_range2));
    waitbar(i / length(valid_alpha_values), h1);
end
close(h1);
% Plot the first radial component (Region 1) after normalization
figure(1);
hold on;
colors = jet(length(valid_alpha_values)); % Color map for eigenvalues
for i = 1:length(valid_alpha_values)
    alpha = valid_alpha_values(i);
    R1_normalized = coefficients1(i) * R1(r_range1, alpha);
    plot(r_range1, R1_normalized, 'Color', colors(i, :), 'LineWidth', 2);
end
xlabel('$r$', 'Interpreter', 'latex', 'FontSize', 16);
ylabel('$R_1(r)$ (Normalized)', 'Interpreter', 'latex', 'FontSize', 16);
title('Normalized Radial Component $R_1(r)$ in Region 1', 'Interpreter', 'latex', 'FontSize', 18);
colorbar;
grid on;
% Plot the second radial component (Region 2) after normalization
figure(2);
hold on;
for i = 1:length(valid_alpha_values)
    alpha = valid_alpha_values(i);
    R2_normalized = coefficients2(i) * R2(r_range2, alpha);
    plot(r_range2, R2_normalized, 'Color', colors(i, :), 'LineWidth', 2);
end
xlabel('$r$', 'Interpreter', 'latex', 'FontSize', 16);
ylabel('$R_2(r)$ (Normalized)', 'Interpreter', 'latex', 'FontSize', 16);
title('Normalized Radial Component $R_2(r)$ in Region 2', 'Interpreter', 'latex', 'FontSize', 18);
colorbar;
grid on;
% Define the full range for plotting
r_range_total = linspace(2 * epsilon, R_star, 200);
% Compute the matched solution R(r)
R_matched_R1 = zeros(size(r_range_total)); % Contribution from R1
R_matched_R2 = zeros(size(r_range_total)); % Contribution from R2
% Loop over each eigenvalue and add contributions from R1 and R2
for i = 1:length(valid_alpha_values)
    alpha = valid_alpha_values(i);
    % Compute the contribution for each range
    R1_contribution = coefficients1(i) * R1(r_range_total(r_range_total <= R0), alpha);
    R2_contribution = coefficients2(i) * R2(r_range_total(r_range_total > R0), alpha);
    % Combine the two contributions
    R_matched_R1(r_range_total <= R0) = R_matched_R1(r_range_total <= R0) + R1_contribution;
    R_matched_R2(r_range_total > R0) = R_matched_R2(r_range_total > R0) + R2_contribution;
end
% Combine the matched solution for the full range
R_matched = R_matched_R1 + R_matched_R2;
% Plot the matched solution
figure(3);
hold on;
% Plot R1 part
plot(r_range_total(r_range_total <= R0), R_matched_R1(r_range_total <= R0), ...
    'LineWidth', 3, 'Color', 'b', 'DisplayName', '$R_1(r)$');
% Plot R2 part
plot(r_range_total(r_range_total > R0), R_matched_R2(r_range_total > R0), ...
    'LineWidth', 3, 'Color', 'r', 'DisplayName', '$R_2(r)$');
% Styling the figure
set(gcf, 'Units', 'Inches', 'Position', [0, 0, 9.8, 7.5], 'PaperUnits', 'Inches', 'PaperSize', [9.8, 7.5]);
ax = gca;
ax.FontSize = 32; 
extraInputs2 = {'interpreter', 'latex', 'fontsize', 32};
xlabel('$r$', extraInputs2{:});
ylabel('$R(r)$ (Matched Solution)', extraInputs2{:});
title('Matched Radial Solution $R(r)$', extraInputs2{:});
legend('show', 'Interpreter', 'latex', 'FontSize', 24, 'Location', 'best');
box on;
ax.LineWidth = 3;
grid on;
hold off;
% Compute and plot the matched full solution p(r, t)
time_values = [0:0.02:0.5, 5:5:30]; % Time points for evaluation
p_total = zeros(length(time_values), length(r_range_total)); % Store solutions
% Progress bar for solution computation
h2 = waitbar(0, 'Calculating matched solutions...');
for t_idx = 1:length(time_values)
    t = time_values(t_idx);
    p_t = zeros(size(r_range_total));
    
    for i = 1:length(valid_alpha_values)
        alpha = valid_alpha_values(i);
        % Add contributions from R1 and R2 in their respective ranges
        p_t(r_range_total <= R0) = p_t(r_range_total <= R0) + ...
            coefficients1(i) * R1(r_range_total(r_range_total <= R0), alpha) .* exp(-alpha^2 * t);
        p_t(r_range_total > R0) = p_t(r_range_total > R0) + ...
            coefficients2(i) * R2(r_range_total(r_range_total > R0), alpha) .* exp(-alpha^2 * t);
    end
    p_total(t_idx, :) = p_t; % Store the matched solution for this time
    waitbar(t_idx / length(time_values), h2);
end
close(h2);
% Plot the matched solution over time
figure(4);
hold on;
colors = jet(length(time_values)); % Color map for time points
for t_idx = 1:length(time_values)
    plot(r_range_total, p_total(t_idx, :), 'Color', colors(t_idx, :), 'LineWidth', 2);
end
% Styling the plot
set(gcf, 'Units', 'Inches', 'Position', [0, 0, 9.8, 7.5], 'PaperUnits', 'Inches', 'PaperSize', [9.8, 7.5]);
ax = gca;
ax.FontSize = 32; 
extraInputs2 = {'interpreter', 'latex', 'fontsize', 32};
xlabel('$r$', extraInputs2{:});
ylabel('$p(r, t)$', extraInputs2{:});
title('Matched Solution $p(r, t)$ over Time', extraInputs2{:});
colorbar;
grid on;
box on;
ax.LineWidth = 3;
hold off;

[Torsten]

Doesn't the normalization condition automatically follow from your boundary conditions at r = 2*eps and r = R* ? I think a second condition at r = R_0 is necessary to fix the solution of your piecewise defined PDE. The usual conditions are continuity of p and continuity of the flux dp/dr at r = R_0.

J and Y seem to be Bessel functions, but could you mention what the function M is ?

[Roi Bar-On]

Thanks for your reply! We don't demand the continuity on the derivative since that the two match distributions we're expecting shouldn't follow first derivative continuty. (we see it from the data we analyze in parallel). Sorry about not providing details about the different functions, this was corrected now and I provided also the MATLAB code. M is the hypergeometric function.. I don't see how the normalization condition stem from the BCs..

[Torsten]

In sum, you need two boundary conditions at r = 2*eps and r = R* plus two conditions at r = R_0 (usually continuity of the variable and its flux, but discontinuous conditions are also allowed). I'm not sure whether your normalization condition is able to be a substitute for a flux condition at r = R_0.

What is your final aim of the computation ? Determine p ? Or determine eigenvalues and eigenfunctions ?

"alpha" in your code is undefined - thus it doesn't work as posted. Please check using the "Run" button after submission.

[Roi Bar-On]

For me it runs smoothly. Pay attention that when you use handle (for example in : R2_squared = @(r) (J0(alpha * r) - (J1(alpha * R_star) / Y1(alpha * R_star)) * Y0(alpha * r)).^2;) , you don't need to define alpha yet..... alpha is first calculated in line 146 : alpha = valid_alpha_values(i); .....

If the error occurs at "Compute the integrals" integral_R1=integral(...);, You can try adding alpha in the handle.. Instead of @(r) have @(r,alpha). It should solve the problem eventhough for me it runs smoothly. I am using R2021b. Maybe the compiler here is not adjusted to work with such handles, in my MATLAB it runs great as I said, sorry for the inconvenience!

Basically I aim to calculate all of what you metioned and that's what the code does.. I'm getting these reults, however as I mentioned, with a relatively small jump betwwen the two different R's and p's....

See the figures I'm getting (attached)

[Torsten]

At the latest, values for alpha are needed here:

% Compute the integrals
integral_R1 = integral(@(r) R1_squared(r) .* exp(-kappa / (2 * D) * r.^2), 2 * epsilon, R0);
integral_R2 = integral(@(r) R2_squared(r) .* r, R0, R_star);

but none are supplied.

As said: Try to adjust the code so that using the "Run" arrow in the menu bar gives the expected results.

[Roi Bar-On]

As I said, they are not computed or needed in the line you state.. It's just a definition in the handle function, not calculations are made yet. This is something technical.. I ran it smoothly as I said on my MATLAB, maybe this compiler here is limited..

[Torsten]

% Compute the integrals
integral_R1 = integral(@(r) R1_squared(r) .* exp(-kappa / (2 * D) * r.^2), 2 * epsilon, R0);
integral_R2 = integral(@(r) R2_squared(r) .* r, R0, R_star);

integral_R1 and integral_R2 are values, not function handles. To compute the values, alpha is needed as a numerical value.

If this piece of code runs with your MATLAB version, I would be deeply surprised:

% Parameters
R0 = 1.6;       % Inner radius
R_star = 2.4;   % Outer radius
D = 0.33;       % Diffusion coefficient
epsilon = 0.4;  % Small cutoff
kappa = 1;      % Reaction coefficient
% Special functions (placeholders for the confluent hypergeometric function M)
M = @(a, b, x) hypergeom(a, b, x); % MATLAB has the hypergeom function
J0 = @(x) besselj(0, x);
J1 = @(x) besselj(1, x);
Y0 = @(x) bessely(0, x);
Y1 = @(x) bessely(1, x);
% Define R1(r)^2
R1_squared = @(r) ( ...
    ( (1 - (alpha.^2 * D) / (2 * kappa)) / 2 * M(2 - (alpha.^2 * D) / (2 * kappa), 3, (kappa / (2 * D)) * (2 * epsilon)^2) ...
    ./ ((alpha.^2 * D) / (2 * kappa) * M(1 - (alpha.^2 * D) / (2 * kappa), 2, (kappa / (2 * D)) * (2 * epsilon)^2)) ...
    .* M(-(alpha.^2 * D) / (2 * kappa), 1, (kappa / (2 * D)) * r.^2) ...
    + M(1 - (alpha.^2 * D) / (2 * kappa), 2, (kappa / (2 * D)) * r.^2) ) ).^2;
% Define R2(r)^2
R2_squared = @(r) (J0(alpha * r) - (J1(alpha * R_star) / Y1(alpha * R_star)) * Y0(alpha * r)).^2;
% Compute the integrals
integral_R1 = integral(@(r) R1_squared(r) .* exp(-kappa / (2 * D) * r.^2), 2 * epsilon, R0);
integral_R2 = integral(@(r) R2_squared(r) .* r, R0, R_star);

[Roi Bar-On]

You're correct....

I have editted the code in the post.. Now it seems that it takes too much time to track even the first eigenvalue (before it was quiet fast).

It ran in my MATLAB as I didn't use ''clear all'' before running the program as probably somewhere in the workspace an array of alphas was there. Sorry

[Torsten]

Do you have a reference where such boundary conditions (dp/dr = 0 at r = 2*eps and r = R*, continuity for p at r = R* and an integral condition) are used to fix a solution for p ? I'm still not convinced that your integral condition can substitute a (usual) condition on the flux dp/dr at r = R_0.