% Example2_12b.m code for Example 2.12b
% Square body, four heat fluxes (X22B11Y22B11T1). Temperatures at middle side points
clear all; format short e
% data:
n=50; sigma=2.0; conduct=52.5; volhtcap=3500000;dtime=.6; L=.03;
q=100000*[1 2 2.5 9.]'; alpha=conduct/volhtcap; dt=alpha*dtime/L^2
% Introduce time steps
di=1:1:n;
times=di*dt;
    ii=1:1:25;
    N=ii.*ii;
% n is for number of time steps. tx is for the times corresponding to X
% Pd0 = Phi(0,t), Pd5 = Phi(0.5,t), P1d = Phi(1,t)
for it=1:n
  for itt=1:4
    tx((it-1)*4+itt)=it*dt;
    end
  A=times(it)*pi^2;
P1d(it)=((-1/6)+times(it)-2*sum(exp(-A*N).*(cos(ii*pi)./N))/pi^2)*L/conduct;
Pd5(it)=((-1/24)+times(it)-2*sum(exp(-A*N).*(cos(ii*pi*.5)./N))/pi^2)*L/conduct;
Pd0(it)=((1/3)+times(it)-2*sum(exp(-A*N)./N)/pi^2)*L/conduct; ttt(it)=times(it);
XX(1,:)=[Pd0(it) Pd5(it) P1d(it) Pd5(it)]; XX(2,:)=[Pd5(it) Pd0(it) Pd5(it) P1d(it)];
XX(3,:)=[P1d(it) Pd5(it) Pd0(it) Pd5(it)]; XX(4,:)=[Pd5(it) P1d(it) Pd5(it) Pd0(it)];
if it == 1
    X=XX;
    else
X=[X; XX];
end
end
T=X*q+20;
rankX=rank(X)
svdX=svd(X)
conX=cond(X)
rng(129);
    y=X*q+sigma*randn(4*n,1); qhat=X\y
q_std = diag( inv( X'*X ) );
p = length(qhat);
CI = tinv( 0.975, 4*n-p ) * sqrt( q_std ) * sqrt( (y - X*qhat)'*(y - X*qhat)/(4*n-p) )
% now repeat with different random data
m = 40;
for im = 1:m
    y=X*q+sigma*randn(4*n,1); qhat=X\y;
    fprintf('%5.0f, %5.3f\n',max( abs((qhat - q))), max( abs((qhat - q)))/CI(1) )
end
