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

 找回密码
 立即注册
搜索
热搜: 活动 交友 discuz
查看: 4051|回复: 0
打印 上一主题 下一主题

svd,Tikhonov正则法误差

[复制链接]

125

主题

209

帖子

3070

积分

版主

Rank: 7Rank: 7Rank: 7

积分
3070
跳转到指定楼层
楼主
发表于 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-11-22 18:01 , Processed in 0.080857 second(s), 19 queries .

Powered by Discuz! X3.2 Licensed

© 2001-2013 Comsenz Inc.

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