Trying to understand the implementation of boundary conditions in this MATLAB code

I am following with the textbook: Lloyd N. Trefethen, Spectral Methods in MATLAB, SIAM and trying to understand a very particular MATLAB code. This problem is shown in program 37 and it considers 2nd order wave equation in 2D with periodic boundary conditions in x and Neumann boundary conditions in y:

and the boundary conditions:

So the code is:

 %x variable in [-A,A], Fourier 
 A=3; Nx=50; dx =2*A/Nx; x= -A+dx*(1:Nx)';
 D2x = (pi/A)^2*toeplitz([-1/(3*(dx/A)^2)-1/6 ...
     0.5*(-1).^(2:Nx)./sin((pi*dx/A)*(1:Nx-1)/2).^2]);

 %y variable in [-1,1], chebyshev:
 Ny=15; 
[Dy,y] = cheb(Ny);
D2y=Dy^2;

 BC=-Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny); %Homogenous Neumann BCs for |y|=1 
 %Grid and initial data:
 [xx,yy] = meshgrid(x,y);
 vv=exp(-8*((xx+1.5).^2 + yy.^2));
 dt = 5/(Nx+Ny^2);
 vvold = exp(-8*((xx+dt+1.5).^2+yy.^2));
 %Time stepping by leap frog forula:
  dt = 5/(Nx+Ny^2);clf; plotgap=round(2/dt); dt=2/plotgap;
 for n=0:2*plotgap
     t=n*dt;
     if rem(n+.5,plotgap)<1
         subplot(3,1,n/plotgap+1), mesh(xx,yy,vv), view(-10,60)
         colormap([0 0 0])
         text(-2.5,1,.5,['t=' num2str(t)],'fontsize',18),,grid off, drawnow
     end
     vvnew= 2*vv - vvold + dt^2*(vv*D2x+D2y*vv);
     vvold = vv; vv = vvnew;
     vv([1 Ny+1],:) = BC*vv(2:Ny,:); %Neumann BCs for |y|=1
     %first row and last row in vv = BC x 
     %From Program 33
 end

Now, I have a similar problem where I am trying to reproduce the exact BCs for a 2D Poisson's equation/problem. What I am struggling to understand is this line:

 BC=-Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny); %Homogenous Neumann BCs for |y|=1 

I am not sure what is the author doing here so I can reproduce these exact BCs in my code. My understanding is we want to replace the first/last rows of D2 with values of first/last rows of D and set them to zero? Meaning we want to set the values if the first derivative to zero? Is this correct or am I misunderstanding the syntax. Thanks.

The Cheb function:

 function [ D, x ] = cheb ( N )
  if ( N == 0 )
    D = 0.0;
    x = 1.0;
    return
  end
  x = cos ( pi * ( 0 : N ) / N )';
  c = [ 2.0; ones(N-1,1); 2.0 ] .* (-1.0).^(0:N)';
  X = repmat ( x, 1, N + 1 );
  dX = X - X';
  
%  Set the off diagonal entries.
  D =( c * (1.0 ./ c )' ) ./ ( dX + ( eye ( N + 1 ) ) );
%  Diagonal entries.
  D = D - diag ( sum ( D' ) );
 
  return
 end

[Torsten]

If this code is part of a book, it should be explained therein how homogeneous Neumann boundary conditions are approximated.

[Umar]

Hi Janee,

You asked, Now, I have a similar problem where I am trying to reproduce the exact BCs for a 2D Poisson's equation/problem. What I am struggling to understand is this line:BC=-Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny); %Homogenous Neumann BCs for y=1 I am not sure what is the author doing here so I can reproduce these exact BCs in my code. My understanding is we want to replace the first/last rows of D2 with values of first/last rows of D and set them to zero? Meaning we want to set the values if the first derivative to zero? Is this correct or am I misunderstanding the syntax

You did ask excellent questions and sounds like you are motivated about implementing Neumann boundary conditions to ensure that the flux or gradient of the solution is appropriately handled at the boundaries, reflecting physical constraints or characteristics of the problem being modeled. Let me break down the line in question and provide a detailed explanation to help you grasp its purpose and functionality.

The line of code you mentioned:

BC = -Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny); %Homogeneous Neumann BCs for y=1

Here's a breakdown of what this line is doing:

1. `Dy([1 Ny+1],[1 Ny+1])`: This part refers to selecting a submatrix from the matrix `Dy`. Specifically, it selects rows 1 and Ny+1, and columns 1 to Ny+1.

2. `\`: This symbol denotes matrix division. In this context, it is used to solve a system of linear equations represented by the two submatrices on either side of the division symbol.

3. `Dy([1 Ny+1],2:Ny)`: This part selects rows 1 and Ny+1 from the matrix `Dy` but only columns from 2 to Ny.

The purpose of this line is to apply homogeneous Neumann boundary conditions for y=1 in the context of your 2D Poisson's equation. By solving this linear system using the specified submatrices from `Dy`will ensure that the first and last rows' derivatives are set to zero, effectively enforcing Neumann boundary conditions at those boundaries.

Regarding your understanding, you are on the right track. The goal is indeed to replace the first and last rows of `D` with values from the first and last rows of `Dy` while setting them to zero. This process ensures that the first derivative values at those boundaries are zero, which aligns with the concept of homogeneous Neumann boundary conditions.

If you have any further questions or need additional clarification on any part of this explanation, feel free to ask.

[Janee]

@Torsten The book explains with "examples" especially with Neumann BCs. There is explaination with other types of BCs such as Dirichlet which is easier in implenetation in my opinion, since it deals with setting the actual values of u(x,y) to zero and not the first derivative.

[Umar]

Hi Janee,

In terms of implementation, you mentioned that Dirichlet boundary conditions are easier in your opinion. This is because they directly deal with setting the function values, as opposed to setting the first derivative. From a technical perspective, implementing Neumann boundary conditions can indeed be more complex due to the need to consider derivatives. However, it's important to note that both types of boundary conditions have their own applications and advantages depending on the specific problem being solved.

However, your preference for Dirichlet boundary conditions due to their perceived ease of implementation is valid. However, it's essential to consider the specific requirements and mathematical nature of the problem at hand when choosing between Neumann and Dirichlet boundary conditions.

Now, dealing with Dirichlet boundary conditions in MATLAB, the focus is on setting the actual values of the function ( u(x, y) ) at the boundaries whic is indeed simpler than handling Neumann boundary conditions that involve setting the first derivative of the function.

By the way @Torsten is also a great contributor mentioned in your comments. My intention was to share knowledge, but I do apologize in advance if I caused any interruption. Let us know if you have further comments or concerns, we will be happy to help you.

[Janee]

@Umar

Thanks a lot for your response, it was very helpful. I attempted to apply my understaning from your posts on a simple 1D example. I think the example works with one issue or "Matrix is close to singular or badly scaled" error. So, my simple 1D example looks like:

clearvars; clc; close all;
N=10;
[D,x] = cheb(N); D2 = D^2;
%Neumann BC at |y|=1
D2([1 N+1],:) = D([1 N+1],:);

a=-1; b=1;
u =-x + x.^3/3 ; %
%u=(1/6)*x .*(6*a*b-3*a*x-3*b*x+2* x.^2);

n =ones(size(u));
dndx = D *n; 
 
dudx = D *u;
prod1 = dndx .* dudx;
du2dx = D * dudx;
prod2 = n .*du2dx;
 
invN = (1./n) ;
 
source =  prod1 + prod2; 
uold = ones(size(u)); 
max_iter =500; 
err_max = 1e-6; 
 
 for iterations = 1:max_iter
    phikmax_old = (max(abs(uold)));
    duoldx =D *uold;
    dnudx = dndx .* duoldx;
    ffh = source; 
    RHS = ffh - dnudx; 
    Stilde =invN .* RHS;

   unew =D2\[0;Stilde(2:N);0];
   
    phikmax = (max(abs(unew)));
 
    if phikmax < err_max 
        it_error = err_max /2;
    else
        it_error = abs( phikmax - phikmax_old) / phikmax;
    end
    if it_error < err_max
        break;
    end
    
    uold = unew;
 
 end

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.297102e-19.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.297102e-19.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.297102e-19.

figure
plot(x,unew-max(unew),'--rs',x,(u-max(u))); 
legend('Num solution','Exact solution')

Above example returns correct results regardless of the error, however, if I apply the same method in the book of Neumann BCs for example:

clearvars; clc; close all;
N=10;
[D,x] = cheb(N); D2 = D^2;

%D2([1 N+1],:) = D([1 N+1],:);
BC=-D([1 N+1],[1 N+1])\D([1 N+1],2:N);
a=-1; b=1;
u =-x + x.^3/3 ; 
%u=(1/6)*x .*(6*a*b-3*a*x-3*b*x+2* x.^2);
 
n =ones(size(u));
dndx = D *n; %
 
dudx = D *u;
prod1 = dndx .* dudx;
du2dx = D * dudx;
prod2 = n .*du2dx;
 
invN = (1./n) ;
 
source =  prod1 + prod2;
 
uold = ones(size(u)); 
max_iter =500; 
err_max = 1e-6;
 
 for iterations = 1:max_iter
    phikmax_old = (max(abs(uold)));
    duoldx =D *uold;
    dnudx = dndx .* duoldx;
    ffh = source; 
    RHS = ffh - dnudx; 
    Stilde =invN .* RHS;
 
   unew = D2\Stilde;
   unew([1 N+1],:) = BC*unew(2:N,:);
   
    phikmax = (max(abs(unew)));
 
    if phikmax < err_max 
        it_error = err_max /2;
    else
        it_error = abs( phikmax - phikmax_old) / phikmax;
    end
    if it_error < err_max
        break;
    end
    
    uold = unew;
 
 end

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.379133e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.379133e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.379133e-17.

figure
plot(x,unew-max(unew),'--rs',x,(u-max(u))); 
legend('Num solution','Exact solution')

So, the above example is not working really and I am not sure if I am still not understaning the original example.

The Cheb function

function [ D, x ] = cheb ( N )
  if ( N == 0 )
    D = 0.0;
    x = 1.0;
    return
  end
  x = cos ( pi * ( 0 : N ) / N )';
  c = [ 2.0; ones(N-1,1); 2.0 ] .* (-1.0).^(0:N)';
  X = repmat ( x, 1, N + 1 );
  dX = X - X';
  
%  Set the off diagonal entries.
  D =( c * (1.0 ./ c )' ) ./ ( dX + ( eye ( N + 1 ) ) );
%  Diagonal entries.
  D = D - diag ( sum ( D' ) );
 
  return
 end

[Umar]

Hi Jane,

It sounds like the code you provided is trying to calculate the Chebyshev differentiation matrix (matlab finction) and solves a differential equation using this matrix using matlab script code that you provided. But don’t underestand the reason behind plotting the results. However, in your code, the warning message “Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.297102e-19”is likely related to the line where you solve the linear system using the backslash operator (\):

unew = D2\Stilde;

To address this issue, you can try using a more robust solver, such as the mldivide function (mldivide(D2, Stilde)), which can handle ill-conditioned matrices more effectively. For more information on this function, please refer to

https://www.mathworks.com/help/matlab/ref/mldivide.html

Additionally, you can try increasing the precision of your calculations by adjusting the tolerance level (err_max) or using a higher precision data type. It's also worth noting that the warning message may not necessarily indicate a problem with your code. In some cases, the matrix may be inherently ill-conditioned due to the nature of the problem being solved. However, it's important to ensure that the results are accurate and reliable.

It may be helpful to provide more information about the specific problem you are trying to solve and any additional constraints or requirements. This will allow for a more targeted and accurate solution.I hope this helps! Let me know if you have any further questions.

[Umar]

Hi Janee,

I was reading your earlier comments regarding, “to reproduce the exact BCs for a 2D Poisson's equation/problem”. So, I did some experiments with my code to solve a 2D Poisson's equation with specific boundary conditions using MATLAB.I started with defining the grid and parameters such as the domain length, number of grid points in x and y directions, and grid spacing. Then, create the grid using meshgrid for x and y coordinates. Define the Poisson's equation function f(x, y) as an example.Set the boundary conditions for the left, right, top, and bottom boundaries. Initialize the solution matrix u with zeros and then, implement the boundary conditions by assigning the specified values to the boundary points. Afterwards, the code solved the Poisson's equation numerically using finite differences for the interior points.Finally, I was able to visualize the solution using surf to plot the 2D surface of the solution. Please see attached snippet code.

Please see attached plot.

[Janee]

@Umar

"It's also worth noting that the warning message may not necessarily indicate a problem with your code"

I believe you are right, sense I am able to "apply" the BCs in a different manner to get rid of the error. I think this solves my issue. However, I tried using my simple 1D example to expand it on a 2D example, but I am failing to move from a 1D to 2D Neumann BCs and having an issue with MATALB syntax.

The 2D problem solves for 2D linear Poisson's equation with periodic BCs along x and Neumann BCs along y and I am trying to take what I learned from 1D Poisson's solver with Neumann BCs and apply it here.

2D example similar to above 1D case:

clearvars; clc; close all;
 
Nx = 2; 
Ny = 10;
Lx =3; 
kx  = fftshift(-Nx/2:Nx/2-1);    % wave number vector
dx = Lx/Nx; % Need to calculate dx
 
ksqu = (sin( kx * dx/2)/(dx/2)).^2 ;
kx = sin(kx * dx) / dx;
%-----------
xi_x =  (2*pi)/Lx;
ksqu4inv = ksqu; 

 %helps with error: matrix ill scaled because of 0s
ksqu4inv(abs(ksqu4inv)<1e-14) =1e-6;
xi  = ((0:Nx-1)/Nx)*(2*pi);
x =  xi/xi_x;
Igl = speye(Ny+1);
div_x_act_on_grad_x = -Igl;  

 
ylow = -1; 
yupp =1; 
Ly = (yupp-ylow);
eta_ygl  = 2/Ly;
[D,ygl] = cheb(Ny);
 
ygl = (1/2)*(((yupp-ylow)*ygl) + (yupp+ylow));
D   = D*eta_ygl;
D2  = D*D;
D2([1 Ny+1],:) = D([1 Ny+1],:); %take 1st and last row of D2 and replace it with first and last row of D 
[X,Y]   = meshgrid(x,ygl);
 
%linear Poisson solved iteratively

%ZNy represents the operation of setting the boundary values of y component
%to zero:
ZNy = diag([0 ones(1,Ny-1) 0]); 
div_y_act_on_grad_y = D2*ZNy;
 
u =-Y + Y.^3/3 ;
uh = fft(u,[],2); 
duxk=(kx*1i*xi_x) .*uh; 
du2xk = (kx*1i*xi_x) .*duxk;
duyk = D *uh; 
du2yk = D *duyk;

n = ones(size(u));

invnek = fft(1./n,[],2);
nh = fft(n,[],2);
dnhdxk = (kx*1i*xi_x) .*nh; 
dnhdyk =D * nh;
 
puhnhk = dnhdxk .* duxk; 
pduhdxdnhdxk = dnhdyk .* duyk;
pduhdx2nhk = nh .* du2xk; 
pnhdudx2k = nh .*du2yk; 
NumSourcek =puhnhk + pduhdxdnhdxk + pduhdx2nhk + pnhdudx2k;
 
uold = ones(size(u));
uoldk = fft(uold,[],2);
err_max =1e-8; 
max_iter = 500;
Sourcek = NumSourcek; 
for iterations = 1:max_iter
    phikmax_old = max(max(abs(uoldk)));
    duhdxk = (kx*1i*xi_x) .*uoldk;
    gradNgradUx = dnhdxk .* duhdxk;
    duhdyk = (D) *uoldk ; 
    gradNgradUy = dnhdyk .* duhdyk; 
    
    RHSk = Sourcek - (gradNgradUx + gradNgradUy);    
    Stilde = invnek .* RHSk;
    for m = 1:length(kx)
        L = div_x_act_on_grad_x * (ksqu4inv(m)*xi_x)+ div_y_act_on_grad_y;
        
        unewh(:,m) = L\(Stilde(:,m)); 
    end
        %enforce BCs, instead of [0;Stilde(2:N);0] in 1D   
        unewh=[zeros(1,Nx);  unewh(2:Ny,:); zeros(1,Nx)];
        
    phikmax = max(max(abs(unewh)));
 
    if phikmax < err_max 
        it_error = err_max /2;
    else
        it_error = abs( phikmax - phikmax_old) / phikmax;
 
 
    end
    if it_error < err_max
        break;
    end
    uoldk = unewh;
 
end
 unew = real(ifft(unewh,[],2)); 
DEsol = unew - u;
DE = max(max(abs(DEsol)));
 
 figure
 surf(X, Y, unew-max(max(unew)));
 colorbar;
 title('Numerical solution of \nabla \cdot (n \nabla u) = s in 2D');
 xlabel('x'); ylabel('y'); zlabel('u_{numerical}');

 figure
 surf(X, Y, u-max(max(u)));
 colorbar;
 title('Exact solution of \nabla \cdot (n \nabla u) = s in 2D');
 xlabel('x'); ylabel('y'); zlabel('u_{exact}');

But this does NOT set the values of the first derivative at +/-1 to zero but the values of u(x,y) to zero instead!! why is that?

The Cheb function

 function [ D, x ] = cheb ( N )
  if ( N == 0 )
    D = 0.0;
    x = 1.0;
    return
  end
  x = cos ( pi * ( 0 : N ) / N )';
  c = [ 2.0; ones(N-1,1); 2.0 ] .* (-1.0).^(0:N)';
  X = repmat ( x, 1, N + 1 );
  dX = X - X';
  
%  Set the off diagonal entries.
  D =( c * (1.0 ./ c )' ) ./ ( dX + ( eye ( N + 1 ) ) );
%  Diagonal entries.
  D = D - diag ( sum ( D' ) );
 
  return
 end

[Umar]

Hi Janee,

Try to modify the differentiation matrix calculation within the Cheb function, (Reason: The Cheb function you are using is setting the values of u(x, y) to zero instead). I just made the modification in your Cheb function, here is revised version,

function [D, x] = cheb(N)

    if (N == 0)
        D = 0.0;
        x = 1.0;
        return
    end

    x = cos(pi * (0:N) / N)';
    c = [2.0; ones(N-1, 1); 2.0] .* (-1.0).^(0:N)';
    X = repmat(x, 1, N + 1);
    dX = X - X';

    % Set the off-diagonal entries.
    D = (c * (1.0 ./ c)') ./ (dX + eye(N + 1));

    % Modify diagonal entries for enforcing boundary conditions on first derivative.

    D(1, :) = zeros(1, N + 1);
    D(N + 1, :) = zeros(1, N + 1);

    % Update diagonal entries.

    D = D - diag(sum(D'));

    return
end

By incorporating modification of diagonal entries for enforcing boundary conditions on first derivative into your Cheb function, you can make sure that the values of the first derivative at +/-1 are correctly set to zero as intended. This modification will help you enforce the appropriate boundary conditions for the first derivative in your numerical solution code.Hope, now all your issues have been resolved. Please let me know if you have any further questions.