% Script: TestSteadyu
% Programmer: J. Mueller and T. Shores
% Rev: 5/19/97  % added economization
% Last Rev: 7/25/01  % added general flux and refinements 
% Description:  Routine for solving a BVP by a pseudo Sinc-Galerkin method.
% The 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
% 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) 
% You must run the file sincfcn or have all these fcns available.
% You must define a test problem by COMPILING TestProbx.m
% This is the routine for problems with *unknown* solution

% GLOBAL PARAMETERS

global tndx = 0;
global hfactor=1.0; % should be between 1 and 2, 1 is default

% 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)
global G;
global B;

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

endfunction

%

function retval = q0p(x)
% usage: q0p(x)
% description: this computes the fcn q0'(x)
global G;
global B;

retval = pp(x)*x*G/B/(x+1)+p(x)*G/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 SYSTEM

function [coef,catzero,c,c0,gamnodes] = bvp(d,mm,alpha,bta)
% usage: [coef,catzero,c] = bvp(d,mm,alpha,bta)
% description: this routine solves the system in question,  using the
% functions defined in this file and return the coefficient matrix, rhs
% vector, and the approximate solution  

% inputs 
% d: width of strip
% mm: starting left index for nodes
% alpha: left exponential rate factor
% bta: right exponential rate factor

% local variables
% nn: finishing index for nodes
% h: step size
% m: size of system
% nodes: vector of nodes
% coef: coefficient matrix
% v0,v1: right hand side vectors
% w0,w1: right hand side parts of soln
% c0: proposed solution node vector for ctilde(x)
% errvec:  error vector of c0-ctilde at nodes

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

% 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);
end
% make sure even I^-1 used
if (rem(mm+nn+1,2)>0)
  nn=nn+1;
end;
m = mm + nn + 1;
nodes = arrayf('psi',[-mm:nn]*hx)';
% set up coefficient matrix
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);
pnodes = arrayf('p',nodes);
ppnodes = arrayf('pp',nodes);
fnodes = arrayf('f',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
%coef1 = hx*diag(Dnodes.*gampnodes.*wpnodes./gamnodes./phipnodes) + ...
%I1*diag(Dnodes.*gampnodes.*wnodes./gamnodes) - ...
%diag(Dnodes.*gamnodes.*wpnodes)*I1*diag(1./gamnodes) - ...
%I1*diag(Dnodes.*gamnodes.*phipnodes.*wnodes)*1/hx*I1*diag(1./gamnodes);  
%coef2 = hx*diag(vnodes.*wpnodes ./phipnodes) + I1*diag(vnodes.*wnodes);
%coef3 = hx*diag(lamnodes .*wnodes./phipnodes);
%coef = coef1 - coef2 + coef3;
% an alternate grouping...made no significant difference, but faster
coef = hx*diag(((Dnodes.*gampnodes./gamnodes-vnodes).*wpnodes + ...
lamnodes.*wnodes)./phipnodes) - ...
(diag(Dnodes.*gamnodes.*wpnodes)+I1*diag( ...
Dnodes.*gamnodes.*phipnodes.*wnodes)/hx)*I1*diag(1./gamnodes) + ...
I1*diag(Dnodes.*gampnodes.*wnodes./gamnodes-vnodes.*wnodes);

% This is the discretization of T(qi,w(x)S(i,hx)(phi(x))), i=1 in v1
v0 = Dnodes.*q0pnodes - vnodes.*q0nodes;
Tq0 = hx*(v0.*wpnodes + lamnodes.*wnodes.*q0nodes)./phipnodes + ...
I1*(v0.*wnodes);
v0 = Dnodes.*q1pnodes - vnodes.*q1nodes;
v1 = hx*(v0.*wpnodes + lamnodes.*wnodes.*q1nodes)./phipnodes + ...
I1*(v0.*wnodes);

% Set up the rhs vectors
v0 = hx*wnodes.*fnodes./phipnodes - Tq0;

% solve system and assign values for the function to return
w0 = coef\[v0,v1];
w1 = w0(:,2);
w0 = w0(:,1);
% Find the approx to ctilde:
sv0 = hx*(lamnodes.*(w0	+ q0nodes) - fnodes)./phipnodes;
sv1 = hx*lamnodes.*(w1 - q1nodes)./phipnodes;
catzero = (D(0)*G/B + sum(sv0))/(v(0)+A*D(0)/B + sum(sv1));
c0 = w0 - catzero*w1;
c = c0 + q0nodes + catzero*q1nodes;


endfunction

%

function retval = residerr(d,mm,alpha,bta)
% usage: residerr(d,mm,alpha,bta)
% description: calls bvp(d,mm,alpha,bta), computes residual 
% error vector

global edecay;
global hfactor;

% compute remaining sinc setup parameters
hx = sqrt(hfactor*pi*d/(alpha*mm));
if edecay
  nn = round(log(alpha*mm*hx/bta)/hx);
else
  nn = round(alpha/bta*mm);
end
% make sure even I^-1 used
if (rem(mm+nn+1,2)>0)
  nn=nn+1;
end;
m = mm + nn + 1;
nodes = arrayf('psi',[-mm:nn]*hx)';

csup0nodes = arrayf('csup0',nodes);
cact = [csoln(0);csup0nodes];

[coef,catzero,c] = bvp(d,mm,alpha,bta);
% form the approximation to the residual
ApproxResid = coef*cact;
% compute the difference with the actual residual
retval = ApproxResid - rhs;

endfunction

%

% An experiment (run a testproblem file first):
sincfcn
% set some parameters here 
alpha = 1;
bta = 1;
d = pi/3;
table = zeros(5,7);
% set up the "solution"
for mm = [150]
  % 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);
  end
  % make sure even I^-1 used
  if (rem(mm+nn+1,2)>0)
    nn=nn+1;
  end;
  % solve the problem and compute the error vector
  [coef,slncatzero,c,c0,gamnodes] = bvp(d,mm,alpha,bta);
  slnmm = mm;
  slnnn = nn;
  slnhx = hx;
disp('Cond of coef:')
disp(cond(coef));
nodes = arrayf('psi',[-mm:nn]*hx)';
disp(norm(arrayf('csoln',nodes)-c));
end
% the main loop
kk = 0;
for mm = [8,16,32,64,128]
  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);
  end
  % make sure even I^-1 used
  if (rem(mm+nn+1,2)>0)
    nn=nn+1;
  end;
  m = mm + nn + 1;
  % solve the problem and compute the error vector
  [coef,catzero,c] = bvp(d,mm,alpha,bta);
  nodes = arrayf('psi',[-mm:nn]*hx)';
  table(kk,1) = mm;
  table(kk,2) = m;
  table(kk,3) = hx;
  table(kk,4) = norm(coef,inf);
  table(kk,5) = norm(inv(coef),inf);
  table(kk,6) = abs(catzero-slncatzero);
  table(kk,7) = norm(c- capprox(slnmm,slnnn,slnhx,c0,gamnodes,slncatzero,nodes),inf);
end;
disp('M_x  m  h  norm(M)^in  norm(M^{-1})_inf |e(0) ||e(x_k)||_inf');
disp(table);