Octave/Matlab script for solving parabolic PDE partial derivate equation, specifically heat equation.

output.txt

$ octave --no-gui heatEquation2.m 
Nx =  10
Nt =  5
dx =  0.10000
dt =  0.0050000
c =  1
xmin = 0
xmax =  1
tmin = 0
tmax =  1
r =  0.50000
rr =    2.2204e-16
Vectors x
ans =

   0.00000
   0.10000
   0.20000
   0.30000
   0.40000
   0.50000
   0.60000
   0.70000
   0.80000
   0.90000
   1.00000

Vector t
ans =

   0.000000
   0.005000
   0.010000
   0.015000
   0.020000
   0.025000

Result traspose(u) 
ans =

 Columns 1 through 8:

   0.00000   0.30902   0.58779   0.80902   0.95106   1.00000   0.95106   0.80902
   0.00000   0.29389   0.55902   0.76942   0.90451   0.95106   0.90451   0.76942
   0.00000   0.27951   0.53166   0.73176   0.86024   0.90451   0.86024   0.73176
   0.00000   0.26583   0.50564   0.69595   0.81814   0.86024   0.81814   0.69595
   0.00000   0.25282   0.48089   0.66189   0.77809   0.81814   0.77809   0.66189
   0.00000   0.24044   0.45735   0.62949   0.74001   0.77809   0.74001   0.62949

 Columns 9 through 11:

   0.58779   0.30902   0.00000
   0.55902   0.29389   0.00000
   0.53166   0.27951   0.00000
   0.50564   0.26583   0.00000
   0.48089   0.25282   0.00000
   0.45735   0.24044   0.00000

Result u - Result of PDE explicit method. 
u =

   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.30902   0.29389   0.27951   0.26583   0.25282   0.24044
   0.58779   0.55902   0.53166   0.50564   0.48089   0.45735
   0.80902   0.76942   0.73176   0.69595   0.66189   0.62949
   0.95106   0.90451   0.86024   0.81814   0.77809   0.74001
   1.00000   0.95106   0.90451   0.86024   0.81814   0.77809
   0.95106   0.90451   0.86024   0.81814   0.77809   0.74001
   0.80902   0.76942   0.73176   0.69595   0.66189   0.62949
   0.58779   0.55902   0.53166   0.50564   0.48089   0.45735
   0.30902   0.29389   0.27951   0.26583   0.25282   0.24044
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000

Result uu - Result of PDE anlytical solution u(x, t) = exp(-pi^2 * t) * sin(pi * x) method. 
uu =

   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.30902   0.29414   0.27997   0.26649   0.25366   0.24145
   0.58779   0.55948   0.53254   0.50690   0.48249   0.45926
   0.80902   0.77006   0.73298   0.69769   0.66410   0.63212
   0.95106   0.90526   0.86167   0.82018   0.78069   0.74310
   1.00000   0.95185   0.90602   0.86239   0.82087   0.78134
   0.95106   0.90526   0.86167   0.82018   0.78069   0.74310
   0.80902   0.77006   0.73298   0.69769   0.66410   0.63212
   0.58779   0.55948   0.53254   0.50690   0.48249   0.45926
   0.30902   0.29414   0.27997   0.26649   0.25366   0.24145
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000

Error between analytical solution and computed PDE explicit metod solution
ans =

   0.0000000   0.0000000   0.0000000   0.0000000   0.0000000   0.0000000
   0.0000000   0.0002451   0.0004665   0.0006657   0.0008446   0.0010045
   0.0000000   0.0004663   0.0008873   0.0012663   0.0016065   0.0019106
   0.0000000   0.0006418   0.0012213   0.0017430   0.0022111   0.0026297
   0.0000000   0.0007545   0.0014357   0.0020490   0.0025993   0.0030914
   0.0000000   0.0007933   0.0015096   0.0021544   0.0027331   0.0032505
   0.0000000   0.0007545   0.0014357   0.0020490   0.0025993   0.0030914
   0.0000000   0.0006418   0.0012213   0.0017430   0.0022111   0.0026297
   0.0000000   0.0004663   0.0008873   0.0012663   0.0016065   0.0019106
   0.0000000   0.0002451   0.0004665   0.0006657   0.0008446   0.0010045
   0.0000000   0.0000000   0.0000000   0.0000000   0.0000000   0.0000000

Maximum absolute error
ans =  0.0032505
Plot results

parabolicPDEHeatEquation.m

% Author: Caio Rodrigues 
%
% Objective solve the parabolic partial equations with 
% the following form, similar to  the heat equation as shown 
% below:
%
% Reference: 
%  Numerical Scientific Computing - Ward Cheney - 6ed page 587
%
%
%    ∂u/∂t = c^2 (∂u2/∂x2)
%
%  Data: 
%    c = 1 ; 
%    dx = 0.1    ; Nx = 10 
%    dt = 0.005  ; Ny =  5
%    Analytical solution: 
%      u(x, t) = exp(-pi^2 * t) * sin(pi * x)
%
% Discretization schema: 
% u[i, j] = u(x[i], t[j])
%  
% where 
%  -> x[i] = x0 + i x dx 
%  -> y[i] = y0 + i x dy 
%
% Note: It was tested on Octave.
%
%=================================================%
clear 
clc 

%----------- Useful functions --------------%

% Analytical PDE solution
f_u_analytical = @(x, t)[ exp(-pi^2 * t) * sin(pi * x)];

function u = applyFn2D (xs, ys, fn)
    nx = length(xs);
    ny = length(ys);
    u = zeros(nx, ny);
    for i = 1:nx
        for j = 1:ny
            u(i, j) = fn(xs(i), ys(j));
        end
    end
end


% Number of intervas of x and t 
Nx = 10
Nt = 5
dx = 0.1
dt = 0.005

% Parameter of the heat equation 
c = 1.0 

% Boundary conditions 
xmin = 0.0 
xmax = 1.0 
tmin = 0.0 
tmax = 1.0 

% Intial condition function   -> u(x, t = 0) = f0(x)
f0 = @(x)[ sin(x * pi) ];
% Boundary Condition function -> u(xmin = 0, t) = g1(t) 
g1 = @(t)[ 0 ];  
% Boundray Condition function -> u(xmax = 1, t) = g2(t)
g2 = @(t)[ 0  ];

% Discretized domain 
%--------------------------------
x = xmin:dx:xmax;
t = tmin:dt:(Nt * dt);

% Discretized solution 
u = zeros(Nx + 1, Nt + 1);
% ---- Initial conditions  
% compute u(xmin, t = 0) or u(i, j = 0)
u(:,1) = f0(x);
% ---- Boundary conditions
% compute u(i = 0, j)
u(1,:)  = g1(t);
% compute u(i = Nx, j) = g2(t[j])
u(Nx+1,:) = g2(t);

% Useful coefficients 
% 
r  = c^2 * dt / dx^2 
rr = 1 - 2 * r
% Assertion for sanity checking %
assert(0 < r && r < 1/2, "Error: for the method to be stable r should between 0 <= r <= 1/5 ")

for j = 1:Nt
    for i = 2:Nx
        u(i, j+1) = r * u(i-1, j) + rr * u(i, j) + r * u(i+1,j);
    end
end

disp('Vectors x')
x'
disp('Vector t')
t'

disp('Result traspose(u) ')
u'

disp('Result u - Result of PDE explicit method. ')
u

disp('Result uu - Result of PDE anlytical solution u(x, t) = exp(-pi^2 * t) * sin(pi * x) method. ')
uu = applyFn2D(x, t, f_u_analytical)

disp('Error between analytical solution and computed PDE explicit metod solution')
abs(u - uu)

disp('Maximum absolute error')
max(max(abs(u - uu)))

disp('Plot results')
figure(1)
mesh(x, t, u')
xlabel('x')
ylabel('t')
zlabel('u(t,x)')
title('PDE plot u(x[i], t[j])')

figure(2)                               
mesh(x, t, uu')
xlabel('x')
ylabel('t')
zlabel('u(t,x)')
title('Analytical solution u(x, t) = exp(-pi^2 * t) * sin(pi * x)')