% Script: TestParab
% Programmer: J. Mueller and T. Shores
% Last Rev: 7/31/01  
% Description: Routine for solving the parabolic eqn:
% rho(x)*c_t - (D(x)*c_x)_x + (v(x)*c)_x + lam(x)*c = f(x,t), 0<x<inf.
%  A*c(0) + B*c_x(0) = G(t), t > 0
%  u(infty,t) = 0.
%
% pbvp is a routine for solving this PDE by a pseudo Sinc-Galerkin 
% method in space and a Crank-Nicolson implicit method in time. 
% pbvp incorporates bvp (reason for incorporation rather than calling
% is efficiency), a routine for solving the steady state ode 
% -(D(x)*c')' + (v(x)*c)' + lam(x)*c = f(x), 0 < x < infty. 
%                   A*c(0) + B*c'(0) = G  
%                           u(infty) = 0 
% 
% The routine bvp is a variational approach is based on this form of 
% the problem:  For all test functions u(x)
%  int(D*c'*u'-v*c*u'-lam*c*u,x=0..infty) = int(f*u,x=0..infty) - 
%    D(0)*G/B*u(0) + (v(0)+D(0)*A/B)*u(0)*c(0) 
%
% The C-N factor is mu, with 0 < mu <= 1. The resulting weighted average
% at sucessive time levels transforms into the variational form 
% Tbar(c^{k+1},u) = R^{k+1}(u) + (1/mu -1)*(R^k(u) - Tbar(c^k,u)), 
% where Tbar differs from our T notation only by replacing the lam of T
% by lam + rho/(mu*dt). 
% You must run the file sincfcn or have all these fcns available.
% You must define a test problem by COMPILING TestParx.m

% GLOBAL SINC PARAMETERS
global alpha;
global bta;
global d;
global mm;
global hx;
global nn;
global m;
global tndx;
global hfactor=1.0; % should be between 1 and 2, 1 is default

% GLOBAL TIME STEPPING PARAMETERS
global klim; % iteration limit
global dt; % time step size
global mu; % C-N parameter

% DEFINE THE WEIGHT FUNCTIONS

function retval = gam(x)
% usage: gam(x)
% description: this computes the weight fcn gam(x) for deriv approx's

%mu = 0.0;
%retval = x^(1+mu)/(10+x)^2;
%retval = x*sech(x);
retval = x/(10+x)^2;

endfunction

%

function retval = gamp(x)
% usage: gamp(x)
% description: this computes the fcn gam'(x)

%mu = 0.0;
%retval = x^mu*((10+x)*(1+mu)-2*x)/(10+x)^3;
%retval = (1-x*tanh(x))*sech(x);
retval = (10-x)/(10+x)^3;

endfunction

%

function retval = w(x)
% usage: w(x)
% description:  this is the integration weight
% Warning: If you change this, you must also change wp

retval = gam(x); 

endfunction

%

function retval = wp(x)
% usage: wp(x)
% description:  this is the derivative of the integration weight

retval = gamp(x);

endfunction

% DEFINE THE TRANSFORMATION FUNCTIONS

function retval = q0(x)
% usage: q0(x)
% description: this computes the fcn q0(x) *without* factor G(t)

global B;

retval = p(x)*x/B/(x+1);

endfunction

%

function retval = q0p(x)
% usage: q0p(x)
% description: this computes the fcn q0'(x) *without* factor G(t)

global B;

retval = pp(x)*x/B/(x+1)+p(x)/B/(x+1)^2;

endfunction

%
function retval = q1(x)
% usage: q1(x)
% description: this computes the fcn q1(x)

global A;
global B;

retval = p(x)*(-(A/B+pp(0)-1)*x+1)/(x+1);

endfunction

%

function retval = q1p(x)
% usage: q1p(x)
% description: this computes the fcn q1'(x)

global A;
global B;

retval = pp(x)*(-(A/B+pp(0)-1)*x+1)/(x+1) + p(x)*(-A/B-pp(0))/(x+1)^2;

endfunction

% DEFINE THE INTEGRANDS FOR BILINEAR FORM T

function retval = T0(x)
% usage: T0(x)
% description: this computes the integrand fcn of T(q0,w*Skh(phi(x)))

global hx;
global tndx;

retval = (D(x)*q0p(x) - v(x)*q0(x))* ... 
(wp(x)*Skh(tndx,hx,phi(x)) + w(x)*phip(x)*Skhp(tndx,hx,phi(x))) + ...
lam(x)*q0(x)*Skh(tndx,hx,phi(x));

endfunction

%

function retval = T1(x)
% usage: T1(x)
% description: this computes the integrand fcn of T(q1,Skh)

global hx;
global tndx;

retval = (D(x)*q1p(x) - v(x)*q1(x))* ... 
(wp(x)*Skh(tndx,hx,phi(x)) + w(x)*phip(x)*Skhp(tndx,hx,phi(x))) + ...
lam(x)*q1(x)*Skh(tndx,hx,phi(x));

endfunction


% DEFINE THE SHIFTED UNKNOWN FUNCTION

function retval = csup0(x)
% usage: csup0(x)
% description: this computes the fcn csup0(x)

retval = csoln(x) - q0(x) - csoln(0)*q1(x);

endfunction

% DEFINE THE CONFORMAL MAPPINGS

function retval = xphi(x)
% usage: phi(x)
% description: this computes the fcn phi(x)

retval = log(sinh(x));

endfunction

%

function retval = xphip(x)
% usage: phip(x)
% description: this computes the fcn phi'(x)

retval = 1/tanh(x);

endfunction

%

function retval = xpsi(x)
% usage: psi(x)
% description: this computes the fcn psi(x)

retval = log(exp(x)+sqrt(exp(2*x)+1));

endfunction


function retval = phi(x)
% usage: phi(x)
% description: this computes the fcn phi(x)

retval = log(x);

endfunction

%

function retval = phip(x)
% usage: phip(x)
% description: this computes the fcn phi'(x)

retval = 1/x;

endfunction

%

function retval = psi(x)
% usage: psi(x)
% description: this computes the fcn psi(x)

retval = exp(x);

endfunction

% BUILD AND SOLVE THE TIME DEPENDENT SYSTEM

function [cnodes,cat0old] = pbvp(d,mm,alpha,bta,klim,dt,mu)
% usage: [cnodes,catzero] = pbvp(d,mm,alpha,bta,klim,dt,mu)
% description: this routine solves the system using C-N time stepping, 
% using the functions defined in this file and return the approximate
% solution at the last time step; assume mu is not 0,1.

% inputs:
% d: width of strip
% mm: starting left index for nodes
% alpha: left exponential rate factor
% bta: right exponential rate factor
% klim: limit on time steps 
% dt: time step size
% mu: C-N parameter

% local variables:
% ?nodes: vectors of nodes of ?
% coef: coefficient matrix
% v0,v1: right hand side vectors
% w0,w1: right hand side parts of soln
% c0: proposed solution node vector for csup0
% errvec:  error vector of c0-csup0 at nodes

% global variables:
global A;
global B;
global edecay;
global hfactor;

% do the initialization code 
% compute remaining sinc setup parameters
hx = sqrt(hfactor*pi*d/(alpha*mm));
% if we have exponential decay, make use of the economization 
if (edecay)
  nn = round(log(alpha*mm*hx/bta)/hx);
else
  nn = round(alpha/bta*mm);
endif
% make sure even I^-1 used
if (rem(mm+nn+1,2)>0)
  nn=nn+1;
endif;
m = mm + nn + 1;
% set up static vectors
nodes = arrayf("psi",[-mm:nn]*hx)';
vnodes = arrayf("v",nodes);
Dnodes = arrayf("D",nodes);
phinodes = arrayf("phi",nodes);
phipnodes = arrayf("phip",nodes);
gamnodes = arrayf("gam",nodes);
gampnodes = arrayf("gamp",nodes);
wnodes = arrayf("w",nodes);
wpnodes = arrayf("wp",nodes);
lamnodes = arrayf("lam",nodes);
rhonodes = arrayf("rho",nodes);
lampnodes = lamnodes + rhonodes/(mu*dt);
lammnodes = lamnodes - rhonodes/((1-mu)*dt);
pnodes = arrayf("p",nodes);
ppnodes = arrayf("pp",nodes);
q0nodes = arrayf("q0",nodes);
q0pnodes = arrayf("q0p",nodes);
q1nodes = arrayf("q1",nodes);
q1pnodes = arrayf("q1p",nodes);

% construct I^{1} just once
I1 = Idn(1,m);

% This is the discretization of T(c^0,S(i,hx)(phi(x))) coef matrix
coefp = hx*diag(((Dnodes.*gampnodes./gamnodes-vnodes).*wpnodes)./ ...
phipnodes) - (diag(Dnodes.*gamnodes.*wpnodes)+I1*diag( ...
Dnodes.*gamnodes.*phipnodes.*wnodes)/hx)*I1*diag(1./gamnodes) + ...
I1*diag(Dnodes.*gampnodes.*wnodes./gamnodes-vnodes.*wnodes); 
coefp = coefp + hx*diag(lampnodes.*wnodes./phipnodes); 
coefm = coefp - hx/(dt*mu*(1-mu))*diag(rhonodes.*wnodes./phipnodes);

% This is the discretization of Tp(q0,w(x)S(i,hx)(phi(x)))
% actually,...not quite for Tq0p,Tq0m...off by factor of G(t)
% to be inserted later as needed.
v0 = Dnodes.*q0pnodes - vnodes.*q0nodes;
Tq0p = hx*(v0.*wpnodes + lampnodes.*wnodes.*q0nodes)./phipnodes + ...
I1*(v0.*wnodes);
Tq0m = hx*(v0.*wpnodes + lammnodes.*wnodes.*q0nodes)./phipnodes + ...
I1*(v0.*wnodes);
v0 = Dnodes.*q1pnodes - vnodes.*q1nodes;
v1p = hx*(v0.*wpnodes + lampnodes.*wnodes.*q1nodes)./phipnodes + ...
I1*(v0.*wnodes);
v1m = hx*(v0.*wpnodes + lammnodes.*wnodes.*q1nodes)./phipnodes + ...
I1*(v0.*wnodes);

% variable initialization for entry into main loop
cmu = (1-mu)/mu;
t = 0;
cat0old = g(0);
Gt = 1;
foldnodes = nodes;
fnewnodes = nodes;
cnodes = arrayf("g",nodes);
c0oldnodes = cnodes - G(0)*q0nodes - cat0old*q1nodes;
c0newnodes = nodes;
for i = 1:m
  foldnodes(i) = f(nodes(i),0);
endfor
% the main time itertion loop
for ii = 1:klim
  t = t + dt;
  for i = 1:m
    fnewnodes(i) = f(nodes(i),t);
  endfor
  % set up the rhs vectors
  v0p = hx*(wnodes ./ phipnodes) .* (fnewnodes + cmu*foldnodes) - ...
  cmu*(coefm*c0oldnodes + G(t-dt)*Tq0m + cat0old*v1m) -G(t)*Tq0p;
  % solve system and assign values for new c0
  w2 = coefp\[v0p,v1p];
  w0 = w2(:,1);
  w1 = w2(:,2);
  % Find the approx to c(0,t):
  sv0=lampnodes.*(G(t)*q0nodes+w0)./phipnodes;
  sv1=lampnodes.*(q1nodes-w1)./phipnodes;
  sf=(fnewnodes + cmu*foldnodes - cmu*lammnodes.*cnodes)./phipnodes;
  cat0new = ((G(t)+cmu*G(t-dt))/hx+sum(sf)-sum(sv0))/sum(sv1);
  c0newnodes = w0 - cat0new*w1;
  % clean up for next iteration
  foldnodes = fnewnodes;
  c0oldnodes = c0newnodes;
  cat0old = cat0new;
  cnodes = c0oldnodes + G(t)*q0nodes + cat0old*q1nodes;
endfor;

endfunction

%

% An experiment follows:
% set the globals here 
alpha = 1;
bta = 1;
d = pi/3;
klimold =20;
mu = 0.5;
sincfcn
kk = 0;
for mm = [32]
  % reset values
  klim = klimold;
  dt = 2/klim;
  table = zeros(3,5);
  % compute remaining sinc setup parameters
  kk = kk + 1;
  % compute remaining sinc setup parameters
  hx = sqrt(hfactor*pi*d/(alpha*mm));
  % if we have exponential decay, make use of the economization 
  if (edecay)
    nn = round(log(alpha*mm*hx/bta)/hx);
  else
    nn = round(alpha/bta*mm);
  endif
  % make sure even I^-1 used
  if (rem(mm+nn+1,2)>0)
    nn=nn+1;
  endif;
  m = mm + nn + 1;
  nodes = arrayf("psi",[-mm:nn]*hx)';
  cnodes = zeros(size(nodes));
  catzero = 0;
  csolnnodes = nodes;
  for ii = 1:m
    csolnnodes(ii) = csoln(nodes(ii),klim*dt);
  end
  for jj = 1:2 
    %save last result
    coldnodes = cnodes;
    coldatzero = catzero;
    % solve the problem and compute the error vector
    [cnodes,catzero] = pbvp(d,mm,alpha,bta,klim,dt,mu);
    table(jj,1) = dt;
    table(jj,2) = mm;
    table(jj,3) = m;
    table(jj,4) = abs(catzero-csoln(0,klim*dt));
    table(jj,5) = norm(cnodes-csolnnodes,inf);
    klim = 2*klim;
    dt = dt/2;
  end
  % now extrapolate
  coldnodes = (4*cnodes - coldnodes)/3;
  coldatzero = (4*catzero - coldatzero)/3;
  table(3,1) = 1e-7;
  table(3,2) = mm;
  table(3,3) = m;
  table(3,4) = abs(coldatzero-csoln(0,klim*dt));
  table(3,5) = norm(coldnodes-csolnnodes,inf);
  disp( ...
  '    dt     M_x         m        |e(0)|    ||e(x_k)||_inf');
  table
end