Numerical solution of 1-dimensional heat equation

Question

I need to solve the following heat equation problem

$\frac{∂u}{∂t} =\frac{∂^2u}{∂x^2}, 0\leq x\leq 1,t\geq 0,$

$u|_{t=0}=f(x), 0\leq x\leq 1,$

$\frac{∂u}{∂x}|_{x=0}=1, t\geq 0,$

$u|_{x=1}=1, t\geq 0,$

Using a scheme: (superscript stands for time iterations, subscript - for space iterations)

$\frac{u^{n+\frac{1}{2}}_j-u^{n}_j}{\tau}=\frac{u^{n}_{j+1}-u^{n}_j-u^{n+\frac{1}{2}}_j+u^{n+\frac{1}{2}}_{j-1}}{h^2}, j=1,...,N-1,$

$\frac{u^{n+1}_j-u^{n+\frac{1}{2}}_j}{\tau}=\frac{u^{n+1}_{j+1}-u^{n+1}_j-u^{n+\frac{1}{2}}_j+u^{n+\frac{1}{2}}_{j-1}}{h^2}, j=1,...,N-1.$

Where should I start the computation process? If i try to start at the beginning of the interval, i get: [with known values put in boxes]

$\frac{u^{\frac{1}{2}}_1-\boxed{u^{0}_1}}{\boxed{\tau}}=\frac{\boxed{u^{0}_2}-\boxed{u^{0}_1}-u^{\frac{1}{2}}_1+u^{\frac{1}{2}}_0}{\boxed{h^2}}$, with two unknowns and the condition $\frac{∂u}{∂x}|_{x=0}=1, t\geq 0$ left unused. (I don't know how to apply this condition)

If i try to start at the end of the interval, i get:

$\frac{u^{\frac{1}{2}}_{N-1}-\boxed{u^{0}_{N-1}}}{\boxed{\tau}}=\frac{\boxed{u^{0}_{N}}-\boxed{u^{0}_{N-1}}-u^{\frac{1}{2}}_{N-1}+u^{\frac{1}{2}}_{N-2}}{\boxed{h^2}}$, with two unknowns,

or:

$\frac{\boxed{u^{\frac{1}{2}}_{N}}-\boxed{u^{0}_{N}}}{\boxed{\tau}}=\frac{{u^{0}_{N+1}}-\boxed{u^{0}_N}-\boxed{u^{\frac{1}{2}}_N}+u^{\frac{1}{2}}_{N-1}}{\boxed{h^2}}$, with one unknown and $u^{0}_{N+1}$ lying outside the interval undefined. (How should i define it, if i should?)

Answer 1

At the left endpoint consider the Taylor polynomial for the first $3$ points: $$\begin{align}u_0&=u_0\\ u_1&=u_0+hu_0^{\prime}+\frac12h^2u_0^{\prime\prime}+O(h^3)\\ u_2&=u_0+2hu_0^{\prime}+2h^2u_0^{\prime\prime}+O(h^3)\end{align}$$ Solving for $u_0^{\prime}$, $$\frac{-3u_0+4u_1-u_2}{2h}=u_0^{\prime}+O(h^2)$$ Since you know $u_0^{\prime}=1$ enforcing this boundary condition at each half time step gives you an equation for $u_0$.

At the right endpoint the equation is $$u_N=u_N$$ Looks tautologial, but since you know $u_N=1$ it gives you an equation for $u_N$ at each half time step. Then at each time step you can find $b_0=u_0^{\prime}=1$, $b_N=u_N=1$ and $$b_j=\frac{u_j^n}{\tau}+\frac{u_{j+1}^n-u_j^n}{h^2}$$ For $1\le j\le N-1$. You have to solve a $3\times3$ system for the first $3$ $u_j^{n+\frac12}$: $$\begin{bmatrix}-\frac3{2h}&\frac4{2h}&-\frac1{2h}\\-\frac1{h^2}&\frac1{\tau}+\frac1{h^2}&0\\0&-\frac1{h^2}&\frac1{\tau}+\frac1{h^2}\end{bmatrix}\begin{bmatrix}u_0^{n+\frac12}\\u_1^{n+\frac12}\\u_2^{n+\frac12}\end{bmatrix}=\begin{bmatrix}b_0\\b_1\\b_2\end{bmatrix}$$ But after that the rest of the system is lower triangular so for $3\le j\le N-1$, $$\left(\frac1{\tau}+\frac1{h^2}\right)u_j^{n+\frac12}=\frac{u_{j-1}^{n+\frac12}}{h^2}+b_j$$ And then $u_N^{n+\frac12}=b_N$.

The second half time time is easier. Let $b_0=u_0^{\prime}=1$, $b_N=u_N=1$, and for $1\le j\le N-1$, $$b_j=\frac{u_j^{n+\frac12}}{\tau}+\frac{u_{j-1}^{n+\frac12}-u_j^{n+\frac12}}{h^2}$$ And since the system is already upper triangular you can solve upwards: $u_N^{n+1}=b_N$ and for $N-1\ge j\ge1$, $$\left(\frac1{\tau}+\frac1{h^2}\right)u_j^{n+1}=\frac{u_{j+1}^{n+1}}{h^2}+b_j$$ And finally $$-\frac{3u_0^{n+1}}{2h}=-\frac{4u_1^{n+1}}{2h}+\frac{u_2^{n+1}}{2h}+b_0$$ Matlab code:

% heat2.m

clear all;
close all;
N = 20;
x0 = 0;
xf = 1;
h = (xf-x0)/N;
x = linspace(x0,xf,N+1);
tau = 0.001;
nsteps = 400;
% boundary conditions
up0 = 1;
uf = 1;
u = zeros(N+1,1); % Initial conditions
u = sin(pi*x/(2*xf));
b = zeros(size(u));
t = 0;
% Enforce initial boundary conditions
u(1) = (4*u(2)-u(3)-2*h*up0)/3;
u(end) = uf;
plot(x,u);
hold on;
title('Solution to heat equation');
xlabel('x');
ylabel('u');
legends = [{['t=' num2str(t)]}];
for k = 1:nsteps,
    b(1) = 1;
    for j = 2:N,
        b(j) = (1/tau-1/h^2)*u(j)+(1/h^2)*u(j+1);
    end
    b(N+1) = uf;
    % 3X3 system to start
    A = [-3/(2*h) 4/(2*h) -1/(2*h);
        -1/h^2 (1/tau+1/h^2) 0;
        0 -1/h^2 (1/tau+1/h^2)];
    % Make A(2,1) = 0
    temp = A(2,1)/A(1,1);
    A(2,:) = A(2,:)-A(1,:)*temp;
    b(2) = b(2)-b(1)*temp;
    % Make A(3,2) = 0
    temp = A(3,2)/A(2,2);
    A(3,:) = A(3,:)-A(2,:)*temp;
    b(3) = b(3)-b(2)*temp;
    % Solve upper triangular system for u(1:3)
    u(3) = b(3)/A(3,3);
    u(2) = (b(2)-A(2,3)*u(3))/A(2,2);
    u(1) = (b(1)-A(1,2)*u(2)-A(1,3)*u(3))/A(1,1);
    % Now we can solve downwards
    for j = 4:N,
        u(j) = (b(j)+u(j-1)/h^2)/(1/tau+1/h^2);
    end
    % The right boundary condition
    u(N+1) = b(N+1);
    b(1) = 1;
    for j = 2:N,
        b(j) = (1/tau-1/h^2)*u(j)+(1/h^2)*u(j-1);
    end
    b(N+1) = uf;
    % Since the system is already upper triangular we can solve from the
    % bottom up
    u(N+1) = b(N+1);
    for j = N:-1:2,
        u(j) = (b(j)+u(j+1)/h^2)/(1/tau+1/h^2);
    end
    % The left boundary condition
    u(1) = (b(1)-4*u(2)/(2*h)+u(3)/(2*h))/(-3/(2*h));
    t = t+2*tau;
    if mod(k,50) == 0,
        plot(x,u);
        legends = [legends {['t=' num2str(t)]}];
    end
end
legend(legends,'Location','best');

Output:

I'm sure I got some indices wrong someplace, but hopefully you get the right general idea.