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

 找回密码
 立即注册
搜索
热搜: 活动 交友 discuz
查看: 3386|回复: 0

svd,Tikhonov正则法误差

[复制链接]

123

主题

207

帖子

2992

积分

版主

Rank: 7Rank: 7Rank: 7

积分
2992
发表于 2016-4-6 11:10:54 | 显示全部楼层 |阅读模式
直接法,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 = 1ositive_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
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Archiver|手机版|小黑屋|52matlab技术网站 ( 粤ICP备14005920号-5 )

GMT+8, 2024-3-29 09:40 , Processed in 0.062524 second(s), 20 queries .

Powered by Discuz! X3.2 Licensed

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表