Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

$$ \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,:)']
share|improve this question

2 Answers 2

up vote 1 down vote accepted

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));
share|improve this answer
    
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
share|improve this answer
    
My bad, I edited the answer. One of indexes was off by one. –  Stuart Mar 4 '13 at 20:01

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.