# Numerical solution of 1-dimensional heat equation

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

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.