5
$\begingroup$

$$ \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,:)']
$\endgroup$
4
$\begingroup$

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));
$\endgroup$
  • $\begingroup$ Thank you @Stuart, now it is ok. $\endgroup$ – Ömer Mar 4 '13 at 20:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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