pm_opposite.m

来自「数值分析, 同济大学教材<<现代数值数学和计算>>数值分」· M 代码 · 共 27 行

function [x, lam, it] = pm_opposite(A, x0, tol, p)
[n m] = size(A);
if n ~= m
  error('Invalid Metrix.');
end
A = A - p * eye(n);
[L U] = lu(A);
it = 0;
v0 = x0;
u = v0 / maxnorm(v0);
y = L \ u;
v = U \ y;
while norm(v - v0) > tol
  it = it + 1;
  if it > 20000
    disp('iteration failed.');
    return;
  end
  v0 = v;
  u = v / maxnorm(v);
  y = L \ u;
  v = U \ y;
end
x = u;
lam = p + 1 / maxnorm(v);