52matlab技术网站，matlab教程，matlab安装教程，matlab下载

标题: svd,Tikhonov正则法误差 [打印本页]

作者: matlab的旋律 时间: 2016-4-6 11:10
标题: svd,Tikhonov正则法误差
直接法，svd和Tikhonov正则法误差分析

clear all
close all
clc
for n = [10,20,50]
[A,b] = Calc_A_b(n);

t = ([1:n]' - 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 [A,b] = 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 = [10,20,50]
[A,b] = Calc_A_b(n);
t = ([1:n]' - 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)
[U,S,V] = svd(A);
Positive_len = length(find(diag(S) > tol));
x = zeros(length(b),1);
for k = 1

ositive_len
x = x + U(:,k)'*b/S(k,k)*V(:,k);
end
end

clear all
close all
clc
miu = 1e-10;%parameter
for n = [10,20,50]
[A,b] = Calc_A_b(n);
t = ([1:n]' - 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)
[U,S,V] = 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

欢迎光临 52matlab技术网站，matlab教程，matlab安装教程，matlab下载 (http://www.52matlab.com/)