# Accuracy of finite difference method for heat equation on a disk

To study an approximation for the heat equation $$\frac{\partial^2 u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\frac{\partial^2 u}{\partial\theta^2}=f(r,\theta)$$ on the disk $D=(0,1)\times(0,2\pi)$ with periodic boundary conditions, we used the following finite difference method $$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta r)^2}+\frac{1}{r_i}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta r}+\frac{1}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta\theta)^2}$$ where $$r_i=(i-\frac{1}{2})\Delta r, \space\space\space\space i=1,2...N+1$$ $$\theta_j=(j-1)\Delta \theta, \space\space\space\space j=1...M+1$$ $$\Delta r=\frac{2}{2N+1}, \space\space\space\space \Delta \theta=\frac{2\pi}{M}$$ $$M=2N$$ By Taylor expansion, we computed the LTE which should be $$\tau=(\Delta r)^2(\frac{1}{12}u^{IV}(r,\theta)+\frac{1}{6}u^{III}(r,\theta))+\frac{(\Delta \theta)^2}{12}u^{IV}(r,\theta)+O((\Delta r)^4,(\Delta \theta)^4)$$ Thus the method should show second order convergence.

We implemented it in MATLAB with the following code (we tried different values of N and different source functions)

N = Nvect(index);

M=2*N;

h = 2/(2*N+1);

k=2*pi/M;

for i=1:N

gv(i) =(i-0.5)*h;

end

for j=1:M

gv2(j)=(j-1)*k;

end

[meshr,meshteta] = meshgrid(gv,gv2);

meshr = meshr';

meshteta = meshteta';

f = @(r,teta) 9.*r;

rhs = f(meshr,meshteta);

u = @(r,teta) r.^3;

u_mesh = u(meshr,meshteta);

boundary=u(gv(N),gv2);

rhs=rhs*(h^2);

F = reshape(rhs,[M*N,1]);

for i=1:N

lambda(i,1)=1/(2*(i-0.5));

end

for j=1:M

F(j*N)=F(j*N)-(1+lambda(N))*boundary;

end

for j=1:N

beta(j,1)=1/((j-0.5)^2*(k^2));

end

I = speye(M);

I2=ones(N,1);

%I2=speye(m);

D=sparse(diag(beta));

d=ones(M);

T=spdiags([1-lambda, -2*I2, 1+lambda],[1,0,-1],N,N)';

Ap= spdiags([d,d,d,d],[1,-1,M-1,-(M-1)],M,M);

A=kron(I,T)-2*kron(I,D)+kron(Ap,D);

U = A\F;

U = reshape(U,N,M);

u_int = u_mesh(1:N,1:M);

e = max(max(abs(U-u_int)));

end


We tried $N = [8,\space16,\space32,\space64,\space128,\space256,\space512,\space1024]$ and the error ratio $\frac{e_{N_{i+1}}}{e_{N_{i}}}$ stays always around 0.5 which shows a first order accuracy. What could cause such inconsistency between theoretical and computational results? What did we get wrong?

EDIT: As a test, we used a source function s.t. it's easy to compute true solution $U$ by hand so we just compute the difference between $U$ and its approximation and take $e_{N_{i}}=max(U_i-U^{appr}_i)$.

• How do you define the error $e_{N_i}$? – nicoguaro Jun 13 '17 at 19:09
• Usually such behaviour is happening because you have one of the indices wrong, which introduces then a first-order error. Btw, do you take into account in error computation that $u_i$ in your case corresponds to $(i-0.5) \Delta r$? – VorKir Jun 14 '17 at 14:40