ecuación del calor con Neumann B. C en matlab

Question

$$ \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 $$ estoy tratando de código de la anterior ecuación del calor con neumann b.c. el uso explícito de diferencias finitas hacia adelante en matlab. He utilizado central de diferencias finitas para las condiciones de contorno. $$ u_{x}(0,t)=\frac{u_{i+1}^{j}-u_{i-1}^{j}}{2h} $$ para i=1 ı utilizado $$ u_{2}^{j}-u_{x}(0,t)2h= u_{0}^{j} $$ y para i=m $$ u_{m-1}^{j}-u_{x}(0,t)2h= u_{m+1}^{j} $$ Aquí mi código y no da los resultados correctos. En realidad no estoy seguro de que me codificado correctamente las condiciones de frontera. Yo llame a la función como heatNeumann(0,0.1,0,1,0.1,0.0001,1) Sería bueno si alguien le puede ayudar.

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,:)']

Answer 1

Usted tiene la idea correcta, su condición de frontera es, $$u_x(t_n,x_0)=\frac{v_1^m-v_{-1}^m}{2h}$$ Ahora a aplicar su plan para obtener $v_0^{m+1}$. Es decir, $$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]$$

Y lo mismo para el derecho de la condición de límite. Luego de su BCs debe ser,

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));

Answer 2

Gracias @Stuart para su responder, pero todavía los resultados numéricos y exacta de los resultados no coinciden. a la izquierda de la columna numérica, columna de la derecha exacta.

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