Why is my code returning incorrect result?

3 views (last 30 days)
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?

Accepted Answer

Aquatris
Aquatris on 16 Feb 2024
Edited: Aquatris on 16 Feb 2024
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
  2 Comments
Janee
Janee on 20 Feb 2024
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
Aquatris on 20 Feb 2024
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!