# heat equation with Neumann B.C in matlab

$$\frac{\partial u}{\partial t}=\alpha\frac{\partial^{2}u}{\partial x^{2}} \qquad u(x,0)=f(x)\qquad u_{x}(0,t)=0\qquad u_{x}(1,t)=2$$ i'm trying to code the above heat equation with neumann b.c. using explicit forward finite differences in matlab. I used central finite differences for boundary conditions. $$u_{x}(0,t)=\frac{u_{i+1}^{j}-u_{i-1}^{j}}{2h}$$ for i=1 ı used $$u_{2}^{j}-u_{x}(0,t)2h= u_{0}^{j}$$ and for i=m $$u_{m-1}^{j}-u_{x}(0,t)2h= u_{m+1}^{j}$$ Here my code and it doesn't give correct results. Actually i am not sure that i coded correctly the boundary conditions. I call the function as heatNeumann(0,0.1,0,1,0.1,0.0001,1) It would be good if someone can help.

function heatNeumann(t_i,t_f,a,b,dx,dt,alpha)

% dt: step size in time dimension
% dx: step size in x axes
% a: left point of domain x
% b: right point of domain x
% alpha: take 1
% t_i: t initial, t_f: t final

n=(t_f-t_i)/dt;
m=(b-a)/dx;
lambda=alpha*dt/dx^2;
T=zeros(m+1,n+1);

x=a:dx:b;
t=t_i:dt:t_f;

u0 = x.^2+1+cos(pi.*x);
T(:,1)=u0; %initial value
u_left = T(2,1); %boundary cond.
u_right = T(m,1)+4*dx; %boundary cond.

for j=1:n
T(1,j+1)=T(1,j)+lambda*(T(2,j)-2*T(1,j)+u_left); %boundary cond.
for i=2:m
T(i,j+1)=T(i,j)+lambda*(T(i+1,j)-2*T(i,j)+T(i-1,j));
end
T(m+1,j+1)=T(m+1,j)+lambda*(u_right-2*T(m+1,j)+T(m,j)); %boundary cond.
u_left = T(2,j+1);
u_right = T(m,j+1);
end

%exact soln

[xx,tt]=meshgrid(x,t);
exact=2.*tt+xx.^2+1+exp(-pi^2.*tt).*cos(pi.*xx);

[T(:,end) exact(end,:)']

-

## 2 Answers

You have the right idea, your boundary condition is, $$u_x(t_n,x_0)=\frac{v_1^m-v_{-1}^m}{2h}$$ Now apply your scheme to get $v_0^{m+1}$. That is, $$v_0^{m+1} = v_0^m + b\mu[v_1^m - 2v_0^m + v_{-1}^m]=v_0^m + b\mu\left[v_1^m - 2v_0^m + \left(v_1^m-2hu_x\left(t_n,x_0\right)\right)\right]$$

And do the same for the right boundary condition. Then your BCs should become,

T(1,j+1)=T(1,j)+lambda*(T(2,j)-2*T(1,j)+T(2,j));
T(m+1,j+1)=T(m+1,j)+lambda*(T(m,j)-2*T(m+1,j)+(T(m,j)+2*dx*2));

-
Thank you @Stuart, now it is ok. –  Ömer Mar 4 '13 at 20:26

Thank you @Stuart for your respond, but still the numerical results and exact results don't match. left column numerical, right column exact.

1.6305    1.5727
1.6123    1.5645
1.5790    1.5415
1.5368    1.5091
1.4946    1.4752
1.4626    1.4500
1.4522    1.4448
1.4745    1.4709
1.5395    1.5385
1.6551    1.6555
1.8264    1.8273

-
My bad, I edited the answer. One of indexes was off by one. –  Stuart Mar 4 '13 at 20:01