VAorthog and VAeval

VAeval.m

@@ -0,0 +1,33 @@
+function [R0,R1] = VAeval(Z,Hes,varargin)  % Vand.+Arnoldi basis construction
+%VAEVAL  Vandermonde with Arnoldi evaluation.
+%
+% This code comes from Brubeck and Trefethen, "Lightning Stokes solver", arXiv 2020.
+% For the mathematics, see Brubeck, Nakatsukasa, and Trefethen, "Vandermonde with
+% Arnoldi", SIAM Review, 2021.
+%
+%  Input:   Z = column vector of sample points
+%         Hes = cell array of Hessenberg matrices
+%         Pol = cell array of vectors of poles, if any
+% Output:  R0 = matrix of basis vectors for functions
+%          R1 = matrix of basis vectors for derivatives
+M = length(Z); Pol = []; if nargin == 3, Pol = varargin{1}; end
+% First construct the polynomial part of the basis
+H = Hes{1}; Hes(1) = []; n = size(H,2); 
+Q = ones(M,1); D = zeros(M,1);
+for k = 1:n
+   hkk = H(k+1,k);
+   Q(:,k+1) = ( Z.*Q(:,k) - Q(:,1:k)*H(1:k,k)          )/hkk;
+   D(:,k+1) = ( Z.*D(:,k) - D(:,1:k)*H(1:k,k) + Q(:,k) )/hkk;
+end
+R0 = Q; R1 = D;
+% Next construct the pole parts of the basis, if any
+while ~isempty(Pol)
+   pol = Pol{1}; Pol(1) = [];
+   H = Hes{1}; Hes(1) = []; np = length(pol); Q = ones(M,1); D = zeros(M,1);
+   for k = 1:np
+      Zpki = 1./(Z-pol(k)); hkk = H(k+1,k);
+      Q(:,k+1) = ( Q(:,k).*Zpki - Q(:,1:k)*H(1:k,k)                   )/hkk;
+      D(:,k+1) = ( D(:,k).*Zpki - D(:,1:k)*H(1:k,k) - Q(:,k).*Zpki.^2 )/hkk;
+   end
+   R0 = [R0 Q(:,2:end)]; R1 = [R1 D(:,2:end)];
+end

VAorthog.m

@@ -0,0 +1,32 @@
+function [Hes,R] = VAorthog(Z,n,varargin)  % Vand.+Arnoldi orthogonalization
+%VAORTHOG  Vandermonde with Arnoldi orthogonalization.
+%
+% This code comes from Brubeck and Trefethen, "Lightning Stokes solver", arXiv 2020.
+% For the mathematics, see Brubeck, Nakatsukasa, and Trefethen, "Vandermonde with
+% Arnoldi", SIAM Review, 2021.
+%
+%  Input:   Z = column vector of sample points
+%           n = degree of polynomial (>= 0)
+%         Pol = cell array of vectors of poles (optional)
+% Output: Hes = cell array of Hessenberg matrices (length 1+length(Pol))
+%           R = matrix of basis vectors
+M = length(Z); Pol = []; if nargin == 3, Pol = varargin{1}; end
+% First orthogonalize the polynomial part
+Q = ones(M,1); H = zeros(n+1,n);
+for k = 1:n       
+   q = Z.*Q(:,k);
+   for j = 1:k, H(j,k) = Q(:,j)'*q/M; q = q - H(j,k)*Q(:,j); end 
+   H(k+1,k) = norm(q)/sqrt(M); Q(:,k+1) = q/H(k+1,k);
+end
+Hes{1} = H; R = Q;
+% Next orthogonalize the pole parts, if any
+while ~isempty(Pol)
+   pol = Pol{1}; Pol(1) = [];
+   np = length(pol); H = zeros(np,np-1); Q = ones(M,1);
+   for k = 1:np       
+      q = Q(:,k)./(Z-pol(k));
+      for j = 1:k, H(j,k) = Q(:,j)'*q/M; q = q - H(j,k)*Q(:,j); end 
+      H(k+1,k) = norm(q)/sqrt(M); Q(:,k+1) = q/H(k+1,k);
+   end
+   Hes{length(Hes)+1} = H; R = [R Q(:,2:end)];
+end

conformal.m

@@ -37,7 +37,9 @@
 % rational functions, Numer. Math. 142 (2019), 359--382.
 %
 % L. N. Trefethen, Numerical conformal mapping with rational 
-% functions, Computational Methods and Function Theory, 2020.
+% functions, Comp. Meth. Funct. Th. 20 (2020), 369-387.
+%
+% See also CONFORMAL2.
 
 %%
 %   Default algorithm: