%  This function solves homework problem [H3.12]
%  ---------------------------------------------
%
echo on
diary nsds_hw3_12.dry
%
%  Start by entering the system
%
A = [ 1250 -25113 -60050 -42647 -23999
       500 -10068 -24057 -17092  -9613
       250  -5060 -12079  -8586  -4826
      -750  15101  36086  25637  14420
       250  -4963 -11896  -8438  -4756 ]
b = [  5 ;  2 ;  1 ; -3 ;  1 ]
c = [ -1 , 26 , 59 , 43 , 23 ]
d = 0
x0 = [ 1 ; -2 ;  3 ; -4 ;  5 ]
%
%  We want to perform three simulations
%
hvec = [ 0.32 , 0.032 , 0.0032 ]
%
%  Let us perform a simulation in Matlab, and use it as
%  the "correct" reference solution, relying on Matlab's
%  capability to solve this problem correctly.
%
h = hvec(3)
S = ss(A,b,c,d)
t = 0:h:10;
u = ones(size(t));
ycorr = lsim(S,u,t,x0);
[n,m] = size(ycorr);
nm = n*m;
%
%  We loop over the three values of h
%
%  Case (a):
%
h1 = hvec(1);
t1 = 0:h1:10;
x = x0;
y1 = c*x0;
[n,m] = size(t1);
n = n*m - 1;
theta = 0.45;
[F,g] = bi45t_lin(A,b,h1,theta);
echo off
for i=1:n,
   x = F*x + g;
   y1 = [ y1 , c*x ];
end
echo on
y = ycorr(1:100:nm)';
%
%  Compute the error
%
errabs = y - y1;
er1 = norm(errabs,'inf')
errrel = errabs ./ y;
erl1 = norm(errrel,'inf')
%
%  Case (b):
%
h2 = hvec(2);
t2 = 0:h2:10;
x = x0;
y2 = c*x0;
[n,m] = size(t2);
n = n*m - 1;
theta = 0.45;
[F,g] = bi45t_lin(A,b,h2,theta);
echo off
for i=1:n,
   x = F*x + g;
   y2 = [ y2 , c*x ];
end
echo on
y = ycorr(1:10:nm)';
%
%  Compute the error
%
errabs = y - y2;
er2 = norm(errabs,'inf')
errrel = errabs ./ y;
erl2 = norm(errrel,'inf')
%
%  Case (c):
%
h3 = hvec(3);
t3 = 0:h3:10;
x = x0;
y3 = c*x0;
[n,m] = size(t3);
n = n*m - 1;
theta = 0.45;
[F,g] = bi45t_lin(A,b,h3,theta);
echo off
for i=1:n,
   x = F*x + g;
   y3 = [ y3 , c*x ];
end
echo on
y = ycorr';
%
%  Compute the error
%
errabs = y - y3;
er3 = norm(errabs,'inf')
errrel = errabs ./ y;
erl3 = norm(errrel,'inf')
%
%  Plot the results
%
subplot(3,1,1)
plot(t,ycorr,'k-',t1,y1,'r-.')
grid on
title('BI4/5(0.45) Linear Integration with Different Step Sizes')
xlabel('Time')
ylabel('y')
subplot(3,1,2)
plot(t,ycorr,'k-',t2,y2,'r-.')
grid on
xlabel('Time')
ylabel('y')
subplot(3,1,3)
plot(t,ycorr,'k-',t3,y3,'r-.')
grid on
xlabel('Time')
ylabel('y')
print -deps nsds_hw3_12.eps
%
%  Close diary and exit
%
diary off
echo off
return