matlab的旋律 发表于 2016-4-6 11:10:54

svd,Tikhonov正则法误差

直接法,svd和Tikhonov正则法误差分析

clear all
close all
clc
for n =
    = Calc_A_b(n);
   
    t = (' - 1)*pi/(n - 1);
    x = A\b;
    figure
    hold on
    plot(t,x)
    plot(0:0.01:pi,sin(0:0.01:pi))
    title('Dirctly approach')
    xlabel(['n = ',num2str(n)])
    hold off
end

function = Calc_A_b(n)
s = ((1:n)' - 1)*pi/(2*(n - 1));
b = (exp(s) - exp(-s))./s;
b(1) = 2;
A = zeros(n,n);
for r = 1:n
    for c = 1:n
      if c == 1 || c == n
            A(r,c) = pi/(2*(n-1))*exp(s(r)*cos((c-1)*pi/(n-1)));
      else
            A(r,c) = 2*pi/(2*(n-1))*exp(s(r)*cos((c-1)*pi/(n-1)));
      end
    end
end
end


clear all
close all
clc
tol = 1e-6;%tolerance
for n =
    = Calc_A_b(n);
    t = (' - 1)*pi/(n - 1);
   x = Calc_x(A,b,tol);
   
    figure
    hold on
    plot(t,x)
    plot(0:0.01:pi,sin(0:0.01:pi))
    title('Truncated SVD approach')
    xlabel(['n = ',num2str(n)])
    hold off
end

function x = Calc_x(A,b,tol)   
= svd(A);
Positive_len = length(find(diag(S) > tol));
x = zeros(length(b),1);
for k = 1:Positive_len
    x = x + U(:,k)'*b/S(k,k)*V(:,k);
end
end


clear all
close all
clc
miu = 1e-10;%parameter
for n =
    = Calc_A_b(n);
    t = (' - 1)*pi/(n - 1);
   x = Calc_x(A,b,miu);

    figure
    hold on
    plot(t,x)
    plot(0:0.01:pi,sin(0:0.01:pi))
    title('Tikhonov regularization approach')
    xlabel(['n = ',num2str(n)])
    hold off
end

function x = Calc_x(A,b,miu)   
= svd(A);
n = length(b);
x = zeros(n,1);
for k = 1:n
    x = x + S(k,k)*U(:,k)'*b/(S(k,k)*S(k,k)+miu)*V(:,k);
end
end
页: [1]
查看完整版本: svd,Tikhonov正则法误差