HeatConduction

% simple model of 1-dimensional heat conduction
%
% Differential equation is:
%   du/dt = k * d^2(u)/dx^2 + f
% where f(x,t) is a heat source (or sink)
%
% We simulate this by using a grid with the second order difference
% approximation:
%  u''(x) = 1/h^2 * (u(x+h) - 2 u(x) + u(x-h)) + O(h^2)

function HeatConduction

k = 1.14;  % thermal conductivity
L = 1;  % length of the rod
nCells = 100;  % number of cells to simulate

leftBoundary = 300;
rightBoundary = 300;

Q = 0.5 * speye(nCells);  % process noise

measuredPoints = 0:.1:1;  % which distances on rod are measured
R = 0.01 * eye(length(measuredPoints));  % measurement noise

initState = [leftBoundary; zeros(nCells-2,1); rightBoundary]; % initial state
%initState = 300 * ones(nCells,1);

initPred = initState;  % initial prediction for state
initCov = zeros(nCells,nCells);  % initial covariance

endTime = 1e-2;

%%%%
%
% model description 

dx = L/(nCells-1);
dt = 1/20 * dx^2/k;  % time increment for discrete simulation
                     % must be less than 1/2 dx^2/k for convergence
                     % since we use an unstable method

diagA(:,1) = 1/dx^2 * ones(nCells,1);
diagA(:,2) = -2/dx^2 * ones(nCells,1);
diagA(:,3) = 1/dx^2 * ones(nCells,1);
A = spdiags(diagA, [-1, 0, 1], nCells, nCells);
F = speye(nCells) + A * dt;
F(1,1) = 1;
F(1,2) = 0;     % fix up boundary conditions to not change
F(end, end-1) = 0;
F(end, end) = 1;

nState = length(measuredPoints);
measuredCells = unique(round(1 + (nCells-1) * measuredPoints));
measuredCells = measuredCells(1 <= measuredCells <= nCells);
H = sparse(1:nState, measuredCells, 1, nState, nCells); % measurement

linSysModel = LinearSystemModel(F, Q);
linObsModel = LinearObserverModel(H, R);
ekfSysModel = EKFSystemModel(linSysModel);
ekfObsModel = EKFObserverModel(linObsModel);
ukfSysModel = UnscentedSystemModel(linSysModel);
ukfObsModel = UnscentedObserverModel(linObsModel);

linFilt = BLUEFilter(linSysModel, linObsModel);
ekfFilt = BLUEFilter(ekfSysModel, ekfObsModel);
ukfFilt = BLUEFilter(ukfSysModel, ukfObsModel);


iter = ceil(endTime/dt) + 1;
stateVec = zeros([nCells, iter]);
stateVec(:,1) = initState;
measVec = zeros([nState, iter]);

trueState = initState;

for k = 2:iter
    % simulate time evolution and measurement
    trueState = linSysModel.simulate(trueState);
    observation = linObsModel.simulate(trueState);
    
    stateVec(:,k) = trueState;
    measVec(:,k) = observation;
end

function predVec = doFilter(filt)
    predVec = zeros([nCells, iter]);
    predVec(:,1) = initPred;
    P = initCov;

    for k = 2:iter
        [predVec(:,k), P] = filt.filter(predVec(:,k-1), ...
                                        P, measVec(:,k));
    end
end

linPredVec = doFilter(linFilt);
ekfPredVec = doFilter(ekfFilt);
ukfPredVec = doFilter(ukfFilt);

times = dt * (0:(iter-1));

figure
colormap('hot');
image(times, [0 L], stateVec, 'CDataMapping', 'Scaled');
title('Dispersion of Heat');
colorbar;

figure
colormap('hsv');
image(times, [0 L], stateVec - linPredVec, 'CDataMapping', 'Scaled');
title('Linear Prediction Model Error');
colorbar;

figure
colormap('hsv');
image(times, [0 L], stateVec - ekfPredVec, 'CDataMapping', 'Scaled');
title('EKF Prediction Model Error');
colorbar;

figure
colormap('hsv');
image(times, [0 L], stateVec - ukfPredVec, 'CDataMapping', 'Scaled');
title('UKF Prediction Model Error');
colorbar;

end