您的位置：首页 > 编程语言 > MATLAB

数值分析: 病态问题 & 算法稳定

2016-11-26 13:36 405 查看

数值分析实验1

1.1 病态问题

实验目的

算法有”优劣”之分, 问题也有”好坏”之别. 在数值分析的问题中, 如线性代数方程组、矩阵特征值、非线性方程及方程组都存在病态问题.

病态问题: 指输出结果对输入数据非常敏感, 如果输入数据有微小误差, 则引起解的误差会非常大.

问题提出

考虑一个高次代数多项式:

p(x)=(x−1)(x−2)...(x−19)(x−20)=∏k=120(x−k) (1)显然, 这个多项式有单重根1,2,3,...20.

现考虑一个扰动:

p(x)+Ex19=0 (2)其中E是非常小的一个数.

分析方程(1)和方程(2)根的差别, 从而得到方程(1)的解对扰动的敏感性.

实验要求

选择充分小的E反复进行上述实验, 记录结果的变化并进行分析.

将扰动项改成Ex18或者其他, 实验中又会有怎么样的现象出现?

Matlab程序

%数值实验1.1病态问题
%输入:[0,20]之间的的扰动项及小的扰动常数
%输出:加扰动后得到的全部根
result = inputdlg({'请选择被扰动项:[0,20]之间的整数:'},'charpt1_1',1,{'19'});
Numb = str2num(char(result));
if((Numb > 20) || (Numb < 0))
errordlg('请选择正确的扰动项:[0,20]之间的整数:');
return;
end
result = inputdlg({'请输入(0,1)之间的扰动常数:'},'charpt1_1',1,{'0.00001'});
ess = str2num(char(result)); %#ok<*ST2NM>
ve = zeros(1,21);
ve(21 - Numb) = ess;
root = roots(poly(1:20)+ve);
disp(['对第',num2str(Numb),'加扰动',num2str(ess),',得到的全部根为:']);
disp(num2str(root));

实验结果

对若干项加扰动0.00001得到的全部根与真实根如下表所示:

高次项有扰动, 对解的影响非常大.( 如何增加算法的抗扰动性?)


扰动第19项	扰动第15项	扰动第10项	扰动第5项	扰动第1项	真实解
22.5961+2.3083i	19.9979	20.0003	20.003	20.0003	20
22.5961-2.3083i	19.0178	18.9972	18.9972	18.9972	19
18.8972+5.00563i	17.919	18.0112	18.0112	18.0112	18
18.8972-5.00563i	17.1646	16.9711	16.9711	16.9711	17
14.9123+4.95848i	15.5656	16.0483	16.0483	16.0483	16
14.9123-4.95848i	15.4595	14.9354	14.9354	14.9354	15
12.0289+3.73551i	13.7373	14.0653	14.0653	14.0653	14
12.0289-3.73551i	13.1828	12.9491	12.9491	12.9491	13
10.059+2.33021i	11.94	12.0334	12.0334	12.0334	12
10.059-2.33021i	11.018	10.984	10.984	10.984	11
8.63829+1.05641i	9.99768	10.0061	10.0061	10.0061	10
8.63829-1.05641i	8.99968	8.99839	8.99839	8.99839	9
7.70894+0i	8.00025	8.00028	8.00028	8.00028	8
7.02801+0i	6.99994	6.99997	6.99997	6.99997	7
5.99942+0i	6.00001	6	6	6	6
5.00001+0i	5	5	5	5	5
4+0i	4	4	4	4	4
3+0i	3	3	3	3	3
2+0i	2	2	2	2	2
1+0i	1	1	1	1	1

1.2 误差传播与算法稳定性

实验目的

算法稳定性非常重要, 误差扩张的算法是不稳定的, 误差衰减的算法是稳定的, 衰减的误差是算法设计的目标.

问题提出

考虑一个简单的定积分序列:

En=∫10xnex−1dx, n=1,2,3...

显然积分序列为递推关系:

En=1−nEn−1, n=1,2,3...

实验要求

我们由递推关系, 可以得到两种算法:

算法1

E1=1e, En=1−nEn−1, n=1,2,3...

算法2

En−1=1−Enn, n=N−1,N−2,...3,2.

分别利用算法1和算法2, 分别采用5、6、7位有效数字, 哪一个能给出精确、稳定的结果.

Matlab程序

%数值实验1.2:误差传播与算法稳定性
%输入:递推公式选择与递推步数
%输出:个布递推值及误差结果,以及递推值和误差与递推步数的关系图
promps = {'请选择递推关系式,若选择算法1,请输入1,否则请输入2'};
result = inputdlg(promps,'charpt1_2',1,{'1'});
Nb = str2num(char(result)); %#ok<*ST2NM>
if((Nb ~= 1) && (Nb ~= 2))
errordlg('请选择递推关系式,若选择算法1,请输入1,否则请输入2!');
return;
end
result = inputdlg({'请输入递推步数n:'},'charpt1_2',1,{'10'});
steps = str2num(char(result));
if(steps < 1)
errordlg('递推步数错误!');
return;
end
result = inputdlg({'请输入计算中所采用的有效数字位数:'},'charpt 1_2',1,{'5'});
Sd = str2num(char(result));
format long %设置显示精度
result = zeros(1,steps);
err = result; %#ok<*NASGU>
func = result;

%用库函数quadl计算积分近似值
for n = 1:steps
fun = @(x) x.^n.*exp(x-1);
func(n) = quadl(fun,0,1); %#ok<*DQUADL>
end

if(Nb == 1)
%用算法1计算
digits(Sd);%控制有效数字位数
result(1) = subs(vpa(1/exp(1)));
for n = 2:steps
result(n) = subs(vpa(1-n * result(n-1)));
end
err = abs(result - func);
elseif(Nb == 2)
%用算法2计算
digits(Sd);%控制有效数字位数
result(steps) = 0;
for n = steps:-1:2
result(n-1) = abs(vpa((1 - result(n)) /n ));
end
err = abs(result - func);
end

clf;%清除当前图像窗口
disp('递推值:');
disp(sprintf('%e  ',result)); %#ok<*DSPS>
disp('误差:');
disp(sprintf('%e  ',err));
plot([1:steps], result,'-'); %#ok<*NBRAK
text(4,result(1),'\downarrow En');
grid on;
hold on;
plot([1:steps], err,'r--');
xlabel('n');
ylabel('En and ERR n');
text(2,err(2),'\uparrow err(n)');

实验结果1

比较算法1的精度, 算法1随着精度次数增加, 误差不稳定.

算法1-递归10-精度5