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,


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


1 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;
hold on;
title('Solution to heat equation');
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);
    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);
    % 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);
    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);
    % 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,
        legends = [legends {['t=' num2str(t)]}];

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


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .