Why is my code returning incorrect result?

I have two codes that solve the same problem, but one is returning the correct solution and the other is not. I don't seem to see where the issue is right away.

Code 1: Solves a 1D linear BVP u_xx = exp(4x), with zero Dirichlet BCs u(a)=u(b)=0

N=10;
[D,x] = cheb(N); 
ylow = 0; %a
yupp = 3; %b
Ly = yupp-ylow;
eta_ygl  = 2/Ly;
x = (1/2)*((yupp-ylow)*x + (yupp+ylow));
D=D*eta_ygl;
D2 = D^2;
D2 = D2(2:N,2:N);
%Exact solution
u = (yupp*exp(4*ylow)-ylow*exp(4*yupp)+ylow*exp(4*x)-yupp*exp(4*x)-exp(4*ylow)*x+exp(4*yupp)*x)/(16*ylow-16*yupp);
n = x.^2+10; 
dndx = D *n; 
dudx = D *u;
prod1 = dndx .* dudx;
du2dx = D * dudx;
prod2 = n .*du2dx;
invN = (1./n) ;
source = prod1(2:N) + prod2(2:N);
uold = ones(size(u(2:N)));
max_iter =500; 
err_max = 1e-4; 

 for iterations = 1:max_iter
    OldSolMax = (max(abs(uold)));
    duoldx = D(2:N,2:N) *uold;
     dnudx = dndx(2:N) .* duoldx;
    RHS = source - dnudx; 
    Stilde = invN(2:N) .* RHS; 

    unew = D2\Stilde;
    NewSolMax = (max(abs(unew)));

    if NewSolMax < err_max 
        it_error = err_max /2;
    else
        it_error = abs( NewSolMax - OldSolMax) / NewSolMax;
    end
    if it_error < err_max
        break;
    end
    uold = unew;
 end
unew = [0;unew;0];
figure
plot(x,unew,'--rs');
legend('Num solution')

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

Code 2: Solves the same problem in Code 1 but it is set in 2D, the code assumes no variations in x and therefore it solves a 1D problem in a 2D setting. This code is essentially the exact same as code 1 and should return the same result:

Nx = 2; 
Ny = 10; 
Lx =3; 
kx  = fftshift(-Nx/2:Nx/2-1);    % wave number vector
dx = Lx/Nx; % Need to calculate dx
% Use approximations for kx, ky, and k^2. These come from Birdsall and Langdon. Find page number and put it here at some point,
ksqu = (sin( kx * dx/2)/(dx/2)).^2 ;
kx = sin(kx * dx) / dx;
xi_x =  (2*pi)/Lx;
ksqu4inv = ksqu; 
ksqu4inv(abs(ksqu4inv)<1e-14) =1; %helps with error: matrix ill scaled because of 0s
xi  = ((0:Nx-1)/Nx)*(2*pi);
x =  xi/xi_x;
ylow = 0;%a
yupp =3;%b 
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;
[X,Y]   = meshgrid(x,ygl);
Igl = speye(Ny+1);
div_x_act_on_grad_x = -Igl;  
%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;
A = 2*pi / Lx;
u = (yupp*exp(4*ylow)-ylow*exp(4*yupp)+ylow*exp(4*Y)-yupp*exp(4*Y)-exp(4*ylow).*Y+exp(4*yupp).*Y)/(16*ylow-16*yupp);
uh = fft(u,[],2);
n = Y.^2+10;
invnek = fft(1./n,[],2);
nh = fft(n,[],2);
dnhdxk = (kx*1i*xi_x).*nh;
dnhdyk =D * nh;
a = ylow; b= yupp;
ExactSource1D =((-exp(4*a) + exp(4*b) + 4*(a - b)*exp(4*Y)) .*Y)/(8*(a - b)) +  exp(4*Y) .*(10 + Y.^2);

uold = ones(size(u));
uoldk = fft(uold,[],2);

err_max =1e-4;
max_iter = 500;
Sourcek = fft(ExactSource1D,[],2); 
for iterations = 1:max_iter
    OldSolMax = max(max(abs(uoldk)));
    duhdxk = (kx*1i*xi_x) .*uoldk;
    gradNgradUx = aapx(dnhdxk,duhdxk);
    duhdyk = (D) *uoldk ; 
    gradNgradUy = aapx(dnhdyk,duhdyk); 
 
    RHSk = Sourcek - (gradNgradUx + gradNgradUy); 
    Stilde = aapx(invnek,RHSk);
    for m = 1:length(kx)
        L = (-Igl* (ksqu4inv(m))*(xi_x^2))+ div_y_act_on_grad_y;
        unewh(:,m) = L\(Stilde(:,m));

    end
      %enforce BCs    
    unewh=[zeros(1,Nx);  unewh(2:Ny,:); zeros(1,Nx)];
    NewSolMax = max(max(abs(unewh)));

    if NewSolMax < err_max 
        it_error = err_max /2;
    else
        it_error = abs( NewSolMax - OldSolMax) / NewSolMax;
    end
    if it_error < err_max
        break;
    end
    uoldk = unewh;

end
 unew = real(ifft(unewh,[],2)); 
 figure
 surf(X, Y, unew);
 colorbar;
 title('Numerical solution of \nabla \cdot (n \nabla u) = s in 2D');

 figure
 surf(X, Y, u);
 colorbar;
 title('Exact solution of \nabla \cdot (n \nabla u) = s in 2D');

aapx function is:

function ph=aapx(uh,vh) %anti-aliased product, 
%removal of aliasing by pading or truncatiion page 134 canto
%use discrete transform with m rather than n points, where m >= 3N/2
[ny nx]=size(uh);m=nx*3/2;
uhp=[uh(:,1:nx/2) zeros(ny,(m-nx)) uh(:,nx/2+1:nx)]; % pad uhat with zeros
vhp=[vh(:,1:nx/2) zeros(ny,(m-nx)) vh(:,nx/2+1:nx)]; % pad vhat with zeros
up=ifft(uhp,[],2); 
vp=ifft(vhp,[],2); 
w=up.*vp; 
wh=fft(w,[],2);
ph=1.5*[wh(:,1:nx/2) wh(:,m-nx/2+1:m)]; % extract F-coefficients 

end 

and Cheb function:

function [ D, x ] = cheb ( N )
%% cheb() computes the Chebyshev differentiation matrix and grid.

  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';
  D = ( c * (1.0 ./ c )' ) ./ ( dX + ( eye ( N + 1 ) ) );
  D = D - diag ( sum ( D' ) );
  return
end

As you can see code 2 is giving a numerical result that is off! The only difference between code 1 and code 2 is the method which the boundray conditions were applied.

Code 1 directly applies BCs on the derivative D,

D2 = D2(2:N,2:N);

where code 2 uses the operator ZNy to applies the BCs on the y compnent of the Laplacian L and re-enforces BCs inside the for loop

ZNy = diag([0 ones(1,Ny-1) 0]);
div_y_act_on_grad_y = D2*ZNy;

and inside the iterative solver:

unewh=[zeros(1,Nx);  unewh(2:Ny,:); zeros(1,Nx)];

How can I fix this so that the two codes give the same results?

Answer 1

[Aquatris]

You might wanna provide the equations yu are trying to implement. Otherwise, unless someone is already familiar with the equations, it is hard to judge what might be wrong. It might be a simple syntax issue (you meant to do a element wise multiplication but did a matrix multiplication) or it can be a more fundemantal issue (you were supposed to select your base from -1 to 1 but you did -1 to 0).

Having said that, i think you issue is here

ksqu4inv(abs(ksqu4inv)<1e-14) =1; %helps with error: matrix ill scaled because of 0s

You do not want to convert 0s to 1s instead you should set 0s to low values like:

ksqu4inv(abs(ksqu4inv)<1e-14) =1e-3; %helps with error: matrix ill scaled because of 0s

to prevent division by 0s

[Janee]

Thanks! That does solve the problem, however, temporarily only! As I increase Ly to 16 for example the solution deviates again and so on! That's why I keep going back to the BCs.

Nx = 2; 
Ny = 10; 
Lx =3; 
kx  = fftshift(-Nx/2:Nx/2-1);    % wave number vector
dx = Lx/Nx; % Need to calculate dx
% Use approximations for kx, ky, and k^2. These come from Birdsall and Langdon. Find page number and put it here at some point,
ksqu = (sin( kx * dx/2)/(dx/2)).^2 ;
kx = sin(kx * dx) / dx;
xi_x =  (2*pi)/Lx;
ksqu4inv = ksqu; 
ksqu4inv(abs(ksqu4inv)<1e-14) =1e-3; %helps with error: matrix ill scaled because of 0s
xi  = ((0:Nx-1)/Nx)*(2*pi);
x =  xi/xi_x;
ylow = 0;%a
yupp =16;%b 
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;
[X,Y]   = meshgrid(x,ygl);
Igl = speye(Ny+1);
div_x_act_on_grad_x = -Igl;  
%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;
A = 2*pi / Lx;
u = (yupp*exp(4*ylow)-ylow*exp(4*yupp)+ylow*exp(4*Y)-yupp*exp(4*Y)-exp(4*ylow).*Y+exp(4*yupp).*Y)/(16*ylow-16*yupp);
uh = fft(u,[],2);
n = Y.^2+10;
invnek = fft(1./n,[],2);
nh = fft(n,[],2);
dnhdxk = (kx*1i*xi_x).*nh;
dnhdyk =D * nh;
a = ylow; b= yupp;
ExactSource1D =((-exp(4*a) + exp(4*b) + 4*(a - b)*exp(4*Y)) .*Y)/(8*(a - b)) +  exp(4*Y) .*(10 + Y.^2);
uold = ones(size(u));
uoldk = fft(uold,[],2);
err_max =1e-4;
max_iter = 500;
Sourcek = fft(ExactSource1D,[],2); 
for iterations = 1:max_iter
    OldSolMax = max(max(abs(uoldk)));
    duhdxk = (kx*1i*xi_x) .*uoldk;
    gradNgradUx = aapx(dnhdxk,duhdxk);
    duhdyk = (D) *uoldk ; 
    gradNgradUy = aapx(dnhdyk,duhdyk); 
 
    RHSk = Sourcek - (gradNgradUx + gradNgradUy); 
    Stilde = aapx(invnek,RHSk);
    for m = 1:length(kx)
        L = (-Igl* (ksqu4inv(m))*(xi_x^2))+ div_y_act_on_grad_y;
        unewh(:,m) = L\(Stilde(:,m));
    end
      %enforce BCs    
    unewh=[zeros(1,Nx);  unewh(2:Ny,:); zeros(1,Nx)];
    NewSolMax = max(max(abs(unewh)));
    if NewSolMax < err_max 
        it_error = err_max /2;
    else
        it_error = abs( NewSolMax - OldSolMax) / NewSolMax;
    end
    if it_error < err_max
        break;
    end
    uoldk = unewh;
end
 unew = real(ifft(unewh,[],2)); 
 figure
 surf(X, Y, unew);
 colorbar;
 title('Numerical solution of \nabla \cdot (n \nabla u) = s in 2D');

aapx

 function ph=aapx(uh,vh) %anti-aliased product, 
%removal of aliasing by pading or truncatiion page 134 canto
%use discrete transform with m rather than n points, where m >= 3N/2
[ny nx]=size(uh);m=nx*3/2;
uhp=[uh(:,1:nx/2) zeros(ny,(m-nx)) uh(:,nx/2+1:nx)]; % pad uhat with zeros
vhp=[vh(:,1:nx/2) zeros(ny,(m-nx)) vh(:,nx/2+1:nx)]; % pad vhat with zeros
up=ifft(uhp,[],2); 
vp=ifft(vhp,[],2); 
w=up.*vp; 
wh=fft(w,[],2);
ph=1.5*[wh(:,1:nx/2) wh(:,m-nx/2+1:m)]; % extract F-coefficients 
 end 
 

Cheb

 function [ D, x ] = cheb ( N )
%% cheb() computes the Chebyshev differentiation matrix and grid.
  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';
  D = ( c * (1.0 ./ c )' ) ./ ( dX + ( eye ( N + 1 ) ) );
  D = D - diag ( sum ( D' ) );
  return
end

[Aquatris]

I think you are lacking understanding on the basic principle of the method and I recommend first understand what each parameter represent and their effects on your solution.

From a first look, I guess that Ly parameter defines something like a stepsize. Becuase you increased it, you lost some info and got a large yellow part in the graph. Check the data points you have in the graph between what you think is wrong and what you think is the right plot and see if they have the same data points but showow different results. I think if you decrease the Ly, than that large yellow part will become more colorful, indicating more data is present there.