Regarding grid generation using finite difference method (*code fixed)

Hi everyone,

I'm trying to generate a uniform grid where the red lines intersect at a 45-degree angle to the horizontal lines in blue using the finite difference method in MATLAB.

However, the resulting grid is slightly off, as shown in the attached image. Intesection angle is uniform but higher than 45-degree.

I set computational domain as and a grid is generated on physical domain with coordinates .

As I want to generate a mesh described above, I aimed to solve the following PDE :

where we have (horizontal lines). The boundary conditions are .

The above equation comes down to the following:

Based on the above equation, I discretized to finite difference form (central difference for and forward differece for ).

Plus, I found the grid collapses when I set different angles (for example 90-degree (replace )).

Am I missing something in implementation?

I need your help.

Here is my code.

% Parameter settings
eta_max = 10;  
xi_max = 10;  
Delta_eta = 0.1;  
Delta_xi = 0.1;  
tol = 1e-6;  
max_iterations = 100000;  

% Grid size
eta_steps = floor(eta_max / Delta_eta) + 1;  
xi_steps = floor(xi_max / Delta_xi) + 1;  

% Initial conditions
x = zeros(eta_steps, xi_steps);  
y = repmat((0:eta_steps-1)' * Delta_eta, 1, xi_steps);  

% Boundary conditions
x(1, :) = linspace(0, xi_max, xi_steps);  

% Finite difference method
converged = false;  
iteration = 0;  

while ~converged && iteration < max_iterations  
    iteration = iteration + 1;  
    max_error = 0;  

    for n = 1:eta_steps - 1  
        for m = 2:xi_steps - 1  
            x_xi = (x(n, m+1) - x(n, m-1)) / (2 * Delta_xi);  
            y_xi = (y(n, m+1) - y(n, m-1)) / (2 * Delta_xi); 

            new_x = x(n, m) + Delta_eta * ((pi/4) - y_xi)/(x_xi);  
            error_x = abs(new_x - x(n, m));  
            max_error = max([max_error, error_x]);  

            x(n + 1, m) = new_x;  
        end  

        x(n+1, xi_steps) = x(n+1, xi_steps-1) + Delta_xi;  
        x(n+1, 1) = x(n+1, 2) - Delta_xi;  
    end  

    if max_error < tol  
        converged = true;  
    end  
end  

% Plot the results
figure; hold on; axis equal; grid on;  

for m = 1:xi_steps  
    plot(x(:, m), y(:, m), 'r');  
end  

for n = 1:eta_steps  
    plot(x(n, :), y(n, :), 'b');  
end  

title(['Iterations until convergence: ', num2str(iteration)]);  
xlabel('x'); ylabel('y');  

Answer 1

[David Goodmanson]

Hello arxiv,

It appears that what you want in the code is not

new_x = x(n, m) + Delta_eta * ((pi/4) - y_xi)/(x_xi);

but rather

new_x = x(n, m) + Delta_eta * (tan(pi/4) - y_xi)/(x_xi);
      = x(n, m) + Delta_eta * (     1    - y_xi)/(x_xi);

which leads to

[ArxIv]

Hi, David Goodmanson.

I replaced the value with and modifed the code a bit.

Finally, I could get the nice grid as attached.

Thank you very much for your comment.