数值实验题

大学的 数值实验题1

实验1.1 病态问题

实验目的：

算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题是指问题本身对扰动敏感，反之属于好问题。本实验通过对一个高次多项式方程的求解，初步认识病态问题。

实验内容：

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

p(x)?(x?1)(x?2)?(x?20)??(x?k) (E.1.1)

k?120显然该多项式的全部根为1，2，?，20，共计20个，且每个根都是单重的(也称为简单的)。现考虑该多项式的一个扰动 p(x)??x19?0， (E.1.2)

其中，ε是一个非常小的数。这相当于是对方程(E.1.1)中x19的系数作一个小的扰动。比较方程(E.1.1)和方程(E.1.2)根的差别，从而分析方程(E.1.1)的解对扰动的敏感性。 实验步骤与结果分析：

(一) 实验源程序

function t_charpt1_1

% 数值实验1.1病态问题

% 输入：[0 20]之间的扰动项及小的扰动常数 % 输出：加扰动后得到的全部根 clc

result=inputdlg({'请输入扰动项:在[0 20]之间的整数:'},'charpt 1_1',1,{'19'}); Numb=str2num(char(result));

if((Numb>20)|(Numb<0))errordlg('请输入正确的扰动项:[0 20]之间的整数!');return;end

result=inputdlg({'请输入(0 1)之间的扰动常数:'},'charpt 1_1',1,{'0.00001'}); ess=str2num(char(result)); ve=zeros(1,21); ve(21-Numb)=ess;

root=roots(poly(1:20)+ve);

x0=real(root); y0=imag(root); plot(x0',y0', '*');

disp(['对扰动项 ',num2str(Numb),'加扰动',num2str(ess),'得到的全部根为:']); disp(num2str(root));

1

(二)实验结果分析

(1) 对于x19项的扰动ess，不同的取值对应的结果如下所示。 ? 对扰动项 19加扰动1e-010得到的全部根为: 19.9961，19.0257，17.9085，17.1508，15.7982，15.181，13.8995，13.0571，11.9753，11.0109，9.99608，9.00111，7.99978，7.00003，6，5，4，3，2，1。 ? 对扰动项 19加扰动1e-009得到的全部根为:

19.952，19.2293，17.6573+0.692896i，17.6573-0.692896i，15.4524+0.875524i，15.4524-0.875524i，13.3527+0.486992i，13.3527-0.486992i，11.8578，11.0427，9.9916，9.00201，7.99952，7.00009，5.99999，5，4，3，2，1。 ? 对扰动项 19加扰动1e-007得到的全部根为:

20.422+0.999203i，20.422-0.999203i，18.1572+2.4702i，18.1572-2.4702i，15.3149+2.69865i，15.3149-2.69865i，12.8466+2.06246i，12.8466-2.06246i，10.9216+1.10366i，10.9216-1.10366i，9.56629，9.11508，7.99387，7.00027，6，5，4，3，2，1。 ? 对扰动项 19加扰动1e-005得到的全部根为:

22.5961+2.3083i，22.5961-2.3083i，18.8972+5.00563i，18.8972-5.00563i，14.9123+4.95848i，14.9123-4.95848i，12.0289+3.73551i，12.0289-3.73551i，10.059+2.33021i，10.059-2.33021i，8.63828+1.0564i，8.63828-1.0564i，7.70896，7.028，5.99942，5.00001，4，3，2，1。

根在复平面上的位置如图所示：

图 ess=1e-010 图 ess=1e-009

图 ess= 1e-007 图 ess=1e-005

2

从实验的图形中可以看出，当ess充分小时，方程E.1.1和方程E.1.2的解相差很小，当ess逐渐增大时，方程的解就出现了病态解，这些解都呈现复共轭性质。并且，病态解首先出现在x=16这个解附近，如ess=1e-009时，x=20，19，12，11，?，2，1的解基本误差不大。在x=16附近，扰动后的解偏离实轴程度较严重，随着ess的增大，扰动对解的影响从x=16附近开始向两边波及，并且偏离实轴的幅度越来越大。x=0，1，2，3，4，5这些阶次较小的解对x19上的扰动最不敏感。

(2)将扰动项加到x18上后，ess=1e-009时方程的解都比较准确，没有出现复共轭现象。 ess=1e-008时误差与x19(ess=1e-009)时相当，即扰动加到x18上比加到x19小一个数量级。 对x8的扰动ess=1000时没有出现复共轭，误差很小；对x的扰动ess=10e10时没有出现复共轭，误差很小。因此，扰动作用到xn上时，n越小，扰动引起的误差越小。

(3) P(x)?(x?1)(x?2)?(x?20)?x20?210x19???20!?0

令 P(x,?)?P(x)??x19?x20?(?210??)x19?? (E.3.1)

P(x,?)的零点均为ε的函数，分别它们记为xi(?),i?1,2,?,20， 显然有

??0,xi(?)?i,i?1,2,,20。研究 xi(?)关于ε的变化情况，将P(x,?)表示为

P(x,?)?(x?x1(?))(x?x2(?))??(x?x20(?))

?(x?xi(?))j?1,j?i?(x?xj(?)) （E.3.2）

?(x?xj(?)) （E.3.3）

dxi(0)。在（E.3.3）两边令ε?0，d?2020xi(?)关于ε的变化或敏感程度可以用其导数表示，故在（E.3.2）两边关于ε求导： x1920d?dxi(?)????(x?xj(?))?(x?xi(?))?d??j?1,j?id???j?1,j?i为了知道原方程的根是如何受扰动ε的影响，需要知道得到 x19dx(0)??id?19?d(x?j)?(x?i)??d?j?1,j?i??20?(x?xj(?))? (E.3.4)

?j?1,j?i?|??0?20在令x?i，则得 idx(0)??id?j?1,j?i?(i?j)，即

,i?1,2,?,20 （E.3.5）

20dxi(0)??d?i19j?1,j?i?(i?j)20由于 xi(?)?xi(0)?dxi(0)dx(0)(??0)，故 |xi(?)?xi(0)|?i?? （E.3.6） d?d?计算表明，对根1，此导数的绝对值只有8.2?10?18，极其微小；但从第7个根起，此导数的绝对值就从2.5?103开始，最大直到2.4?109，非常大！所以必定造成病态。这就是根源。

现在来估计扰动对根的影响。对根i,i?1,2,?,20的影响，由（E.3.5）可得条件数

3

|?|| Condi(?)?dxi(0)|d? （E.3.7） |i|经过计算发现，扰动对x16?16的影响最大，对x1?1的影响最小,与实际计算结果一致。

实验1.2 误差传播与算法稳定性

实验目的：

体会稳定性在选择算法中的地位。误差扩张的算法是不稳定的，是我们所不期望的；误差衰竭的算法是稳定的，是我们努力寻求的。 实验内容：

考虑一个由积分定义的序列 En??xnex?1dx,n?1,2?, (E.1.4)

01显然En>0，n=1，2，?当n=1时，利用部分积分得E1=1/e。

而对n≥2，经分部积分可得递推关系 En=1-nEn-1, n=2,3,? (E.1.5) 由式(E.1.4)得 En≤1/(n+1)。 由递推关系式(E.1.5)，可得计算式(E.1.4)积分序列{En}的两种算法。其一为式(E.1.5)的直接应用，即

E1=1/e, En= 1 – nEn-1, n=2, 3, ? (E.1.7)

另一种算法则是利用式(E.1.5)变形得到

EN=0，En-1=(1- En)/n, n= N-1, N-2, ?, 3, 2. (E.1.8)

实验步骤及结果分析：

(一)实验源程序

function t_charpt1_2

% 数值实验1.2:误差传播与算法稳定性 % 输入:递推式选择及递推步数

% 输出：各步递推值及误差结果，以及递推值和误差与递推步数的关系图 clc

promps={'请选择递推关系式，若选E.1.7,请输入1，否则输入2:'}; result=inputdlg(promps,'charpt 1_2',1,{'1'}); Nb = str2num(char(result));

if((Nb~=1)&(Nb~=2)) errordlg('请选择递推关系式，若选E.1.7,请输入1，否则输入2!');return;end

result=inputdlg({'请输入递推步数n:'},'charpt 1_2',1,{'10'}); steps=str2num(char(result));

if(steps<1)errordlg('递推步数错误!');return;end

result=zeros(1,steps); err=result; if(Nb==1)

n=1; result(n)=1/exp(1); while(n

4

result(n+1)=1-n*result(n);

err(n+1)=abs(result(n+1)-func(n+1)); n=n+1; end

elseif(Nb==2) n=steps;

err(n)=abs(result(n)-func(n)); while(n>1)

result(n-1)=(1-result(n))/n;

err(n-1)=abs(result(n-1)-func(n-1)); n=n-1; end end

disp(['递推值:',num2str(result)]); disp(['误差:',num2str(err)]); plot([1:steps],result,'-'); grid on hold on;

plot([1:steps],err,'r--');

xlabel('n'); ylabel('En- and ERRn--'); text(2,err(2),'\%uparrow err(n)'); text(4,result(4),'\\downarrow En');

%------------------------------------------------------------------ function en=func(n) % 计算En的精确值 if(n==1)

en=1/exp(1); else

en=1-n*func(n-1); end

(二) 实验结果分析

(1) 分别用算法(E.1.7)、算法(E.1.8)计算得到的结果如图所示。

图 算法(E.1.7)结果 图 算法(E.1.8)结果

5

本文来源：https://www.bwwdw.com/article/ij22.html