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
endI 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 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 neumann MATLAB Answers — New Questions
