I am trying to solve for a 1D BVP with Neumann boundary conditions. I am using chebyshev derivative and so to implement these BCs I am initially replacing the first/last row of the second derivative to the values of first/last row of the first derivative and set them to zero so that u'(-1)=u'(1)=0 (i.e. first derivative at left and right boundries = 0). I am not sure if that’s what I am doing in the code though because I get completely wrong numerical solution.
clearvars; clc; close all;
N=10;%

[D,x] = cheb(N); D2 = D^2;
%D2(1,:) = D(1,:); %Right
%D2(N+1,:) = D(N+1,:); %Neumann BC at left endpoint (x=-1), esp w/ D
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 =(1/6)*x .*(6*a*b-3*a*x-3*b*x+2* x.^2);
%Exact Solution
figure
plot(x,u);
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-8;

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 = D2Stilde; %
%Enforce BCs
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

figure
plot(x,unew,’–rs’,x,(u));
legend(‘Num solution’,’Exact solution’)
Function Cheb is
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 trying to solve for a 1D BVP with Neumann boundary conditions. I am using chebyshev derivative and so to implement these BCs I am initially replacing the first/last row of the second derivative to the values of first/last row of the first derivative and set them to zero so that u'(-1)=u'(1)=0 (i.e. first derivative at left and right boundries = 0). I am not sure if that’s what I am doing in the code though because I get completely wrong numerical solution.
clearvars; clc; close all;
N=10;%

[D,x] = cheb(N); D2 = D^2;
%D2(1,:) = D(1,:); %Right
%D2(N+1,:) = D(N+1,:); %Neumann BC at left endpoint (x=-1), esp w/ D
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 =(1/6)*x .*(6*a*b-3*a*x-3*b*x+2* x.^2);
%Exact Solution
figure
plot(x,u);
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-8;

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 = D2Stilde; %
%Enforce BCs
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

figure
plot(x,unew,’–rs’,x,(u));
legend(‘Num solution’,’Exact solution’)
Function Cheb is
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 trying to solve for a 1D BVP with Neumann boundary conditions. I am using chebyshev derivative and so to implement these BCs I am initially replacing the first/last row of the second derivative to the values of first/last row of the first derivative and set them to zero so that u'(-1)=u'(1)=0 (i.e. first derivative at left and right boundries = 0). I am not sure if that’s what I am doing in the code though because I get completely wrong numerical solution.
clearvars; clc; close all;
N=10;%

[D,x] = cheb(N); D2 = D^2;
%D2(1,:) = D(1,:); %Right
%D2(N+1,:) = D(N+1,:); %Neumann BC at left endpoint (x=-1), esp w/ D
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 =(1/6)*x .*(6*a*b-3*a*x-3*b*x+2* x.^2);
%Exact Solution
figure
plot(x,u);
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-8;

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 = D2Stilde; %
%Enforce BCs
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

figure
plot(x,unew,’–rs’,x,(u));
legend(‘Num solution’,’Exact solution’)
Function Cheb is
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, boundary, condition MATLAB Answers — New Questions

​