Implementing Finite Difference Scheme

I am trying to implement a finite difference scheme where the nodes are h distance apart.

To solve for v

The domain is discretized so that each of the nodes are h distance apart.

However, the boundary conditions lie h/2 away from the last node in each direction.

The boundary condition on the left node is Neumann and Dirichlet on the right node.

I am having troulbe with the left most node, it does not appear to be giving me the correct solution.

The right most node seems to be ok.

Any help would be appreciated.

%% Discretize the Domain
N=10; % Number of nodes in r and z direction
ar=0;
br=1;
dr=1/N;
r=(ar:dr:br-dr)'; % Edge Centered 
rc=.5*(r(1:end-1)+r(2:end)); %Cell Centered 
lambda = (dr./(2*rc)); % needed for differential in R direction on Cell Centered Grid
%Test Case #1
fxn = @(r) 16*r.^2;
uv= @(r) r.^3   ;
%These stay the same no matter what the test function is
rhs = fxn(rc);
vexact = uv(rc);
fv = dr^2*rhs;
% Apply Boundary Conditons
fv(1)=0; % Neumann B.C
fv(end) = fv(end)-((40*lambda(end)+48)/15)*uv(r(end)); % Dirichlet B.C.
m=N-1;
upper_diag_v = 1+lambda(2:end)';
lower_diag_v = 1-lambda(2:end)';
upper_diag_v=[0 0   upper_diag_v]; 
lower_diag_v=[lower_diag_v 0 0];
e=ones(1,m);
Ap1 = spdiags(upper_diag_v'.*e ,1,m,m);
Am1 = spdiags(lower_diag_v'.*e,-1,m,m);
e=ones(m,1);
Ad = spdiags(-2*e,0,m,m);
A=(Am1+Ad+Ap1);
A(1,1)=-2*dr;
A(1,2)=3*dr;%(1+lambda(1));
A(1,3)=-dr;
A(end,end)=-(2*lambda(end)+5);
A(end,end-1)=(6-2*lambda(end))/3;
A(end,end-2)=-1/5;
vapprox = A\fv;

[Torsten]

Please show your discretization using the mathematical formulae. It's tedious to translate your code into the underlying differencing scheme.

Answer 1

[William Rose]

@Kenneth,

I ran your code. It runs without error.

Please describe the problem you are solving in words, and specify the equation you are solving. The script says the nodes extend in the r and z direction, but z is not defined. zMin, zMax, and dz are not defined in the script. That surprises me.

This appears to be a problem involving axisymmetric flow in an annular cylinder with inner diameter ar and outer diameter br. Since ar=0, it is a cylinder without a hole in the middle. The symmetry of the problem means that the Dirichlet boundary condition of zero derivative at the center (at r=0) makes sense for some problems, such as dv/dr=0 at r=0.

You mention problems at the right side and the left side. What do you mean by left and right? Left is z=0, right is z=L?

Please explain the boundary conditions in words.

What are the physical meanings of fv, rhs, vexact, A, fxn(), uv()? Please specify the units for each. This will help others understand the problem, and may help identify inconsistencies.

You have the same number of nodes in the r and z direction. I recommend making the numbers of nodes different, to prevent confusion about which dimension of arrays correspond to each physical direction.

After the script runs, A is a 9x9 array which is mostly tri-diagonal, but it has 4 elements in column 3, and four elements in column 7. Does that suprise you?

After the script runs, lambda=[1, 1/3, 1/5, 1/7, ..., 1/17]'. Is that what you expect?

After the script runs, fv=

[0 0.0036 0.0100 0.0196 0.0324 0.0484 0.0676 0.0900 -2.3316]'.

Is that what you expect?

Answer 2

[Torsten]

Approximate the derivatives in r, not in rc.

Use formula (3.2 c) from the attached document (with c = 1 in your case) for r = 0 and (3.1 c) for inner grid points.