%% Fitting an artificial frequency response % Mario Berljafa \and Stefan Guettel % % February 2015 % % Tags: RKFIT, rational approximation %% Introduction % This is an example of RKFIT being used for approximating an artifically % created frequency response function $f(z)$ and finding its poles. % RKFIT is described in [2, 3] and this code reproduces the example % from [2, Section 5.1]. % The rational function $f(z)$ is taken from [5, Section 4.1]. %% Problem setup % We first define a rational function $f(z)$ via its poles, residues, and % the linear term. The numbers are taken from [5]. pls = [-4500, -41000, -100+5000i, -100-5000i, ... -120+15000i, -120-15000i,-3000+35000i,-3000-35000i, ... -200+45000i, -200-45000i,-1500+45000i,-1500-45000i, ... -500+70000i, -500-70000i,-1000+73000i,-1000-73000i, ... -2000+90000i,-2000-90000i]; rsd = [-3000, -83000, -5+7000i, -5-7000i, ... -20+18000i, -20-18000i, 6000+45000i,6000-45000i, ... 40+60000i, 40-60000i, 90+10000i, 90-10000i, ... 50000+80000i,50000-80000i, 1000+45000i,1000-45000i, ... -5000+92000i,-5000-92000i]; pls = pls(:); rsd = rsd(:); dd = 0.2; hh = 2e-5; f = @(zz) arrayfun(@(z) sum(rsd./(z-pls)) + dd + z*hh,zz); %% % The frequency range of interest is $0$ to $100$ kHz and we can plot the % response function easily. The plot should look exactly the same as % [5, Figure 1], provided the number of sampling points $N$ is reduced % to $100$ (in this case some peaks are missed in the plot). figure(1) N = 200; zz = 1i*linspace(-1e5, 1e5, N).'; plot(abs(zz)/1e4, abs(f(zz)), 'k-') %% Testing RKFIT % We now run RKFIT from [1] to compute a rational least squares approximant to $f(z)$. % To this end we need to "rewrite" the least squares problem in terms of % a matrix function acting on a vector, i.e., $\| F \mathbf b - r(A)\mathbf b \|_2/\|F \mathbf b\|_2 \rightarrow \min$. % The fitting points |zz| will be the eigenvalues of a diagonal matrix $A$, % and $F$ will be a diagonal matrix of the fitting values. The vector $\mathbf b$ can % be used for weighting the data points and we will set all weights equal % to one. A = spdiags(zz, 0, N, N); F = spdiags(f(zz), 0, N, N); b = ones(N, 1); %% % We also need to specify the initial poles of $r(z)$. By choosing $18$ initial % poles the rational function to be computed by RKFIT would be of type % $(18,18)$. However, the above rational function $f(z)$ is of type $(19,18)$, % so we need to increase the numerator degree by $k=1$. This parameter $k$, % together with the maximal number of iterations to perform, % can be given to RKFIT in form of a |param| structure. init = logspace(3,5,9); % initial poles init = [init , conj(init)]; param.k = 1; param.maxit = 10; %% % That's it. Now we can run RKFIT. [xi, ratfun, misfit_rkfit] = rkfit(F, A, b, init, param); %% % Here is a convergence plot of RKFIT, showing the relative least squares % approximation error after each iteration. figure(2) semilogy(misfit_rkfit, 'ro--') hold on xlabel('iteration') title('relative 2-norm misfit') legend('RKFIT') %% The |rkfun| class % The second output |ratfun| is an object of class |rkfun|. As expected % it represents a rational function of type $(19,18)$: disp(ratfun) %% % This rational function $r(z)$ can be evaluated in two ways: % The first option is to evaluate a matrix function $r(A)\mathbf b$ by % calling |ratfun(A,b)| with two input arguments. For example, here we are % calculating the relative least squares fitting error, which will coincide % (up to roundoff errors) with the misfit of the last RKFIT iterate: format shorte disp([norm(F*b - ratfun(A, b))/norm(F*b), misfit_rkfit(end)]) %% % Alternatively, we can evaluate $r(z)$ pointwise by giving only one input % argument. Let's compare $r(z)$ to $f(z)$ plotted before: figure(1), hold on plot(abs(zz)/1e4, abs(ratfun(zz)), 'r--') xlabel('imag(z)/1e4') title('abs(function value)') legend('exact','RKFIT') %% % One may be interested in the poles of $r(z)$ computed by RKFIT. The poles % of an |rkfun| can be easily retrieved using the |poles| method. Let us % compare |poles(ratfun)| to the exact poles of the function $f(z)$. figure(3) plot(pls,'ko'), hold on plot(poles(ratfun),'r*') axis([-5e4, 1e4, -1e5, 1e5]) legend('exact poles','RKFIT poles','Location','NorthWest') %% % The agreement between the poles is very good. %% Some other choices for the initial poles % In practice one may not know the number of poles required to fit % a function to a certain accuracy. % Let us see what happens if we use too many, in this case $24$, initial % poles for RKFIT. init24 = logspace(6,9,12); init24 = [init24 , conj(init24)]; % Rerun RKFIT and plot convergence. [xi, ratfun, misfit_rkfit24] = rkfit(F, A, b, init24, param); figure(2), semilogy(misfit_rkfit24, 'rs--') legend('RKFIT, 18 init poles' , 'RKFIT, 24 init poles'); %% % We may also initialize RKFIT with % $18$ infinite initial poles, which is a convenient choice if no % a priori information about the location of the poles is available. initinf = repmat(inf, 1, 18); % Rerun RKFIT and plot convergence. [xi, ratfun, misfit_rkfitinf] = rkfit(F, A, b, initinf, param); semilogy(misfit_rkfitinf, 'r*--') xlabel('iteration') title('relative 2-norm error') legend('RKFIT, 18 init poles' , ... 'RKFIT, 24 init poles',... 'RKFIT, 18 infinite init poles'); axis([1, 10, 1e-15, 1e2]) %% % For all three choices of initial poles we find that RKFIT converges % within $1$ or $2$ iterations. In fact, Corollary 3.2 in [3] asserts that % RFKIT converges in a single iterations if $F$ is a rational function. % It is due to the presence of rounding errors that the convergence is % delayed slightly. %% Comparison with vector fitting % We now compare RKFIT with the vector fitting code described in % [4, 5]. % % RKT_BIGBREAK % % The following code makes sure the |vectfit3| implementation of % VFIT is on MATLAB's search path, and if it is, defines some % options. if exist('vectfit3') ~= 2 warning('Vector fitting VFIT3 not found in the Matlab path.'); warning('VFIT can be downloaded from:') warning('http://www.sintef.no/Projectweb/VECTFIT/Downloads/VFUT3/') warning('Skipping comparison with VFIT3.'); return end xi_vfit = init; % 18 Initial poles from above opts.relax = 1; % Relaxed non-triviality constraint. opts.stable = 0; % Do not enforce stable poles. opts.asymp = 3; % Options are 1, 2, or 3. opts.skip_pole = 0; opts.skip_res = 0; opts.spy1 = 0; opts.spy2 = 0; opts.logx = 0; opts.logy = 0; opts.errplot = 0; opts.phaseplot = 0; %% % We now run $10$ iteration for the $18$ and $24$ finite poles used % before. The VFIT implementation does not allow for infinite % poles, as those would lead to Vandermonde matrices. % Use the 18 finite poles as initial guess for VFIT. for iter = 1:10 [SER, xi_vfit, rmserr, fit] = ... vectfit3(f(zz.'), zz.', xi_vfit, b.', opts); misfit_vfit(iter) = norm(f(zz.') - fit)/norm(F*b); end ffit = @(zz) arrayfun(@(z) sum((SER.C).'./(z-diag(SER.A))) ... + SER.D + z*SER.E,zz); figure(1) plot(abs(zz)/1e4, abs(ffit(zz)),'b:') xlabel('imag(z)/1e4') title('abs(function value)') legend('exact', 'RKFIT', 'VFIT', 'Location', 'NorthWest') figure(2), hold on semilogy(misfit_vfit, 'bo:') legend('RKFIT', 'VFIT') figure(3) plot(xi_vfit,'bx') legend('exact poles', 'RKFIT poles', 'VFIT poles', ... 'Location', 'NorthWest'); % Use the 24 poles as initial guess for VFIT. xi_vfit = init24; for iter = 1:10 [SER, xi_vfit, rmserr, fit] = ... vectfit3(f(zz.'), zz.', xi_vfit, b.', opts); misfit_vfit24(iter) = norm(f(zz.') - fit)/norm(F*b); end % And plot the errors. figure(2) semilogy(misfit_vfit24,'bs:') legend('RKFIT, 18 init poles', ... 'RKFIT, 24 init poles', ... 'RKFIT, 18 infinite init poles', ... 'VFIT, 18 init poles', ... 'VFIT, 24 init poles'); %% % *Note:* _The above warnings outputted by VFIT are due to ill-conditioned % linear algebra problems being solved, a problem that is circumvented with % RKFIT by the use of discrete-orthogonal rational basis % functions._ % % RKT_BIGBREAK % % This is the end of this example. The following creates a thumbnail. figure(1), plot(NaN) %% References % % [1] M. Berljafa and S. Guettel. _A Rational Krylov Toolbox for % MATLAB,_ MIMS EPrint 2014.56 (), % Manchester Institute for Mathematical Sciences, % The University of Manchester, UK, 2014. % % RKT_BIGBREAK % % [2] M. Berljafa and S. Guettel. _Generalized rational Krylov % decompositions with an application to rational % approximation,_ SIAM J. Matrix Anal. Appl., 36(2):894--916, % 2015. % % RKT_BIGBREAK % % [3] M. Berljafa and S. Guettel. % _The RKFIT algorithm for nonlinear rational approximation,_ % SIAM J. Sci. Comput., 39(5):A2049--A2071, 2017. % % RKT_BIGBREAK % % [4] B. Gustavsen. _Improving the pole relocating properties of % vector fitting,_ IEEE Trans. Power Del., 21(3):1587--1592, % 2006. % % RKT_BIGBREAK % % [5] B. Gustavsen and A. Semlyen. _Rational approximation of % frequency domain responses by vector fitting,_ IEEE % Trans. Power Del., 14(3):1052--1061, 1999.