PLA(Padé近似)技术在二次特征值问题中的应用

PAL for QEP with low-rank damping

一类特殊的二次特征值问题：

• 二次特征值问题（QEP）
  
  Q(λ)x=(λ2M+λC+K)x=0)
• 低秩属性：
  
  rank(C)=r<<n

PAL算法：

• 给定一个位移σ和Padé近似的阶m。
• 计算比例因子。
• 通过某种线性化技术，将问题转化为
  
  Ls(μ)xL=(As−μBs)xL=0
• 计算以上问题的第k小的特征对(μ^,x^L)，使其后向误差满足一定的精度条件。
• 返回QEP问题第k个近似特征对，以及相应的后向误差。

下面给出该算法的代码：

%   Demo routine for PAL methods.%   Testing problem: 200-by-200 damped_beam. %   This is an illustration Examples 1 in%  %clear; close all;% User define parameterssigma = 1e6i;   % shiftm = 1;          % Pade order% Load test data load damped_beam_200.matn = length(K);% Matrix norm computingnm = norm(M,1);nc = norm(C,1);nk = norm(K,1);% Rank revealing decomposition C = E Fl = 1; % rankE = zeros(n,1);E(n/2,1) = sqrt( nc );F = E;% Sigma split sigma = sigma1 * sigma2sigma1 = sqrt(sigma);sigma2 = sqrt(sigma);% Pade approximation parametersd = (2 * m + 1);gamma = 2 / d  * sin([1:m]' * pi / d ).^2;  xi = cos([1:m]' * pi / d).^2;a = (gamma./xi).^.5;D = - diag(xi);poles = - 1 ./ xi;% PAL scaling factorsz1 = abs(sigma)^2 * nm;z2 = 2 * m * abs(sigma) * nc;zeta = 1 / max([z1, z2, nk]);sqt_zeta = sqrt(zeta);% Form the LEP (A, B)M_sigma = - sigma^2 * M;K_sigma = K - M_sigma;E_sigma1 = sigma1 * E * ( kron(eye(l), a') );F_sigma2 = sigma2 * F * ( kron(eye(l), a') );lm = l*m;As = [  zeta * (K_sigma + sigma * d * C),    sqt_zeta * E_sigma1;        sqt_zeta * F_sigma2.',              eye(lm)];Bs = [  zeta * M_sigma,     zeros(n,lm);        zeros(lm,n),        kron(eye(l), D)];nL = length(As);% Compute eigenvalue by eig[V MU] = eig(As,Bs);% Recover eigenvalues, and sort according to distance to sigmamu = diag(MU);Lam_Q = sigma * sqrt(mu + 1);[~, idx] = sort(abs(Lam_Q - sigma), 'ascend');mu = mu(idx);Lam_Q = Lam_Q(idx);V_LEP = V(:,idx);V_Q = V_LEP(1:n, :);  % Compute relative backward errors of the QEPres_PAL = zeros(nL, 1);  for j = 1:nL    e = Lam_Q(j);    v = V_Q(:,j);    v = v / norm(v);    res_PAL(j)  = norm( e^2 * M * v + e * C * v + K * v ) /...        (abs(e)^2 * nm + abs(e) * nc + nk);    V_Q(:, j) = v;end% Compute amplification factors.% idxx = [1,61,154,2,57,153]; % index for 6 eigenvalues of interest idxx = [1,61,152,2,57,151]; % index for 6 eigenvalues of interest (only for m = 1)alphas = zeros(1,6);emus = zeros(1, 6);betas = zeros(1,6);cx = zeros(1,6);alpha_eta_beta = zeros(1, 6);res_LEP = zeros(6,1);nA = norm(As,1);nB = norm(Bs,1);phis = @(e) nA + abs(e) * nB;Gs = @(e) [eye(n), - sqt_zeta * E_sigma1/(eye(lm) - e * kron(eye(l), D))];rho = @(e) nm * abs(e)^2 + nc * abs(e) + nk;mu_e = @(e) sqrt(e + 1) + a'*((eye(m) - e * D )\a) - d;for j = 1:6    jj = idxx(j);    mus = mu(jj);    lams = Lam_Q(jj);    x_Ls = V_LEP(:, jj); nx_Ls = norm(x_Ls);    x = x_Ls(1:n); nx = norm(x);    rhos = rho(lams);    res_LEP(j) = norm(As * x_Ls - mus * Bs * x_Ls) / nx_Ls / (nA + abs(mus) * nB);    alphas(j) = norm(Gs(mus)) * phis(mus) * nx_Ls / nx / zeta / rhos;    emus(j) = abs( mu_e(mus) );    ncx(j) = norm( C * x) / nx;    betas(j) = abs(sigma) * emus(j) * ncx(j) / rhos;    alpha_eta_beta(j) = alphas(j) * res_LEP(j) + betas(j);end% Output resultsfprintf(' Relative residual norm relations\n');fprintf('#  %15s  %15s %8s  %8s  \n', 'Real(lambda)', 'Imag(lambda)',...     'res_Ls', 'res_Q');for j = 1:6     jj = idxx(j);    lams = Lam_Q(jj);    fprintf('%1d    %+1.6e    %1.6e    %1.2e    %1.2e \n',...         j,real(lams),imag(lams), res_LEP(j), res_PAL(jj));end fprintf('\n Error bounds\n');fprintf('#  %8s  %9s  %8s  %8s %8s \n', 'alphas', 'e(mus)',...     '|Cx|/|x|', 'betas', 'alphas*etaLs+betas');for j = 1:6     jj = idxx(j);    lams = Lam_Q(jj);    fprintf('%1d    %4.3f    %3.2e    %3.2e    %3.2e    %3.2e\n',...         j,alphas(j), emus(j), ncx(j), betas(j), alpha_eta_beta(j));end % ================================================================% Eig plot% ================================================================  % Compute eig for Pade order 3% Pade approximation parametersm = 3;d = (2 * m + 1);gamma = 2 / d  * sin([1:m]' * pi / d ).^2;  xi = cos([1:m]' * pi / d).^2;a = (gamma./xi).^.5;D = - diag(xi);poles = - 1 ./ xi;% PAL scaling factorsz1 = abs(sigma)^2 * nm;z2 = 2 * m * abs(sigma) * nc;zeta = 1 / max([z1, z2, nk]);sqt_zeta = sqrt(zeta);% Form the LEP (A, B)M_sigma = - sigma^2 * M;K_sigma = K - M_sigma;E_sigma1 = sigma1 * E * ( kron(eye(l), a') );F_sigma2 = sigma2 * F * ( kron(eye(l), a') );lm = l*m;As = [  zeta * (K_sigma + sigma * d * C),    sqt_zeta * E_sigma1;        sqt_zeta * F_sigma2.',              eye(lm)];Bs = [  zeta * M_sigma,     zeros(n,lm);        zeros(lm,n),        kron(eye(l), D)];% Compute eigenvalue mu by eig, and recover Lambda[V MU] = eig(As,Bs);mu = diag(MU);Lam_Q3 = sigma * sqrt(mu + 1);% Compute 'exact' eig by polyeig + PLV scalingomega = sqrt(nk / nm);zeta = 2/(nk + omega * nc);Lam = polyeig(zeta*K, omega*zeta*C, omega^2*zeta*M);Lam = Lam * omega;figure;subplot(1,2,1)plot(Lam, '.k'); xlim([-15,2]);legend('exact','Location','NorthWest');xlabel('real(\lambda)');ylabel('imag(\lambda)');rectangle('Position',[-8 .6E6 9 1.9E6]);subplot(1,2,2);plot(Lam, '.k'); hold on;plot(Lam_Q, 'ok');plot(Lam_Q3, 'sk');axis([-8, 1, .6E6, 2.5E6]);legend('exact', 'PAL m=1', 'PAL m=3');xlabel('real(\lambda)');ylabel('imag(\lambda)');return;

更细致理论，可以参考：点击进入

参考文献：A Padé approximate linearization algorithm for solving the quadratic eigenvalue problem with low-rank damping