% script: example1.m
% use mu(t)=  (1-exp(-malpha*t))/malpha for this example
% also use F(x,t)=0, alpha(t)=0 (the rhs of right bdy cnd),
% phi(x)=2-x/2+x^2*sin(3*x)/3 and c(x)=1+cos(3*x)/3.
% Check files alpha.m, F.m and mu.m.  Other functions
% are hardwired into this example.
% Choices of example 1 are malpha=0.2 and malpha=3.0 for Example 5.1
% To get the noise and attempted smoothing of Example 5.2 deremark the
% appropriate lines below.

global mindex;
global malpha;
% files work for Octave or Matlab

malpha = 3.0;
el = 2 % right hand spatial limit
T = 2 % time limit
gmma = 1.0 % gamma parameter as above
bta = 1.0 % beta boundary parameter as above
M = 200 % space node indices from dx..M
N =  200 % time node indices from dx..N
dx = el/M;
xnodes = linspace(dx,el,M)';
pnodes = ones(size(xnodes))+.3*cos(3*xnodes);
% here we wire d into qnodes...we're assuming d=1.
qnodes = 2*ones(size(xnodes))-xnodes/2+xnodes.^2.*sin(3*xnodes)/3;
xnodesh = linspace(dx/2,el,2*M)';
pnodesh = ones(size(xnodesh))+.3*cos(3*xnodesh);
qnodesh = 2*ones(size(xnodesh))-xnodesh/2+xnodesh.^2.*sin(3*xnodesh)/3;
qnodesorig = qnodes;

% Here is an operator experiment:
mindex = 0; % now mu calculates mu(t)
dx = el/M;
x = linspace(dx,el,M)';
% this is our "measurement"
[gnodes,Unodes]=forwardsolvrex(M,M,el,T,gmma,bta,pnodes,qnodes,pnodesh,qnodesh);
gnodes = [mu(T);gnodes];
% try randomizing effect for Example 5.2
%gnodes = gnodes+0.0001*(rand(size(gnodes))-0.5);
% try smoothing this effect
%gnodes = polyval(polyfit([0;xnodes],gnodes,10),[0;xnodes]);
%gtmp=[gnodes(1);gnodes;gnodes(M+1)];
%gnodes = 1/3*(gtmp(1:M+1)+gtmp(2:M+2)+gtmp(3:M+3));
% and fix the boundary values
% now construct a first approximation to q(x), namely h(x)
mindex = 1; % now mu calculates mu'(t)
mumax = 0;
tmpt = linspace(0,T,N);
for jj = 1:N
  mumax = max(mumax,abs(mu(tmpt(jj))));
end;
qnodes = ((1+pnodes(1:M-1)*dx/2).*gnodes(1:M-1) - ...
         2*gnodes(2:M) + (1-pnodes(1:M-1)*dx/2).*gnodes(3:M+1));
% this calculation assumes d=1, else stick it in as factor
qnodes = [qnodes;0]/(mumax/gmma*(1-exp(-gmma*T))*dx*dx);
% do the last entry by quadratic interpolation
qnodes(M) = 3*qnodes(M-1)-3*qnodes(M-2)+qnodes(M-3);
plot(xnodes,qnodes);
disp('qnodes(5):')
disp(qnodes(5))
hold on
for jj = 1:8
  qtmp = operatorA(M,N,el,T,gmma,bta,gnodes,pnodes,qnodes);
  norm(qtmp-qnodesorig,inf)
  qnodes = qtmp;
  if (jj < 4)
    plot(x,qnodes);
    disp('qnodes(5):')
    disp(qnodes(5))
  end;
end;
% finally, plot the true value
plot(x,qnodesorig)
disp('qnodesorig(5):')
disp(qnodesorig(5))
% what does this computed q(x) result in? Let's see...
mindex = 0; % mu now calculates mu(t)
newgnodes = [mu(T);forwardsolv(M,N,el,T,gmma,bta,pnodes,qnodes)];
disp('Here is the error in original g minus implied g:');
disp(norm(gnodes-newgnodes,inf));