数值分析上机实验

刘力辉 201021080401

1.“画圆为方”问题也是古希腊人所提出几何三大难题中的另一个问题。即求作一个正方形，使其面积等于已知圆的面积。不妨设已知圆的半径为 R = 1，试用数值试验显示“画圆为方”问题计算过程中的误差。 （1）MATLAB 程序：

y=pi^(1/2); % to generate 15-bit value of square root of pi b=1; d=1; for k=1:8 b=b*10;

d=d/10; % b and d combined to control the digit of x x=d*fix(b*y); s(k)=x^3; l(k)=x; end

format long [l',s']

（2）误差分析：

2． 设，I x x

d x n n

=

+⎰

50

1

由 x n = x n + 5 x n – 1 – 5 x n – 1 可得递推式 I n = – 5I n – 1 + 1/ n

（1）从 I 0 尽可能精确的近似值出发，利用递推公式：

I I n

n n =-+-51

1 ( n = 1,2,…20)

计算从 I 1 到 I 20 的近似值；

（2）从 I 30 较粗糙的估计值出发，用递推公式：

I I n

n n -=-+1151

5 ( n= 30,29, …, 3, 2 )

计算从I 1 到 I 20 的近似值；

(3) 分析所得结果的可靠性以及出现这种现象的原因。 I 0 =dx x

⎰+1

051=ln (5+x )1

0|=ln 6-ln 5 所以I0≈0.18232155679395

format long

I0=log2(6)/log2(exp(1))-log2(5)/log2(exp(1)) % calculate the value of I0=ln6-ln5 for n=1:20

I0=-5*I0+1/n; % recycling equation between I(n+1) and I(n) s(n)=I0; end s'

则计算结果为：

表1

I1 0.0883922160302300 I11 0.0140713362538500 I2 0.0580389198488700 I12 0.0129766520640700 I3 0.0431387340890000 I13 0.0120398166027400 I4 0.0343063295550100 I14 0.0112294884148600 I5 0.0284683522249700 I15 0.0105192245923700 I6 0.0243249055418100 I16 0.0099038770381400 I7 0.0212326151481100 I17 0.0093041442210800 I8 0.0188369242594600 I18 0.0090348344501700 I9 0.0169264898137900 I19 0.0074574066965100 I10 0.0153675509310500

I20 0.0127129665174600

从计算的数据看出I 20=0.0127129665174600 > I 19=0.0074574066965100

又I n 的积分范围为0～1，所以应该有I n >I n+1。所以算法不稳定。 下面分析导致算法不稳定的原因： 令S n 为近似值，则有

n S S n n /151+-=- (1)

又

n I I n n /151+-=- (2)

(1)-(2)得

S n -I n =-5S n-1-I n-1=(-5)n (S 0-I 0)=(-5)n ε0 (3) ε0为I 0的误差。可知由下往上递推时误差是以指数增加传递的，即 εn =(-5)n ε0

所以必然会导致算法不稳定。又公式(3)知，当由上往下递推时，误差是以指数减少传递的，即

ε0=1/(-5)n εn (4) 所以此时算法是稳定的。下面便给出由I 30的估计值算I 29～I 1的值。

⎰

1

6

n

x ≤I x x d x n n =+⎰501≤dx x n ⎰105 所以有

)1(61+n ≤I n ≤)

1(51+n 即1/186≤I 30≤1/155，取为1/155。

format long

I0=1/155 % calculate the value of I30=ln6-ln5

for n=30:-1:2

I0=(-1/5)*I0+1/(5*n); % recycling equation between I(n+1) and I(n) s(n-1)=I0;

end

j=1;

for i=29:-1:1

a(j)=s(i);

j=j+1;

end

a'

计算出相应数据为：

表2

S29 0.005376344086020 S14 0.011229186626480 S28 0.005821282906930 S13 0.012039876960420 S27 0.005978600561470 S12 0.012976639992530 S26 0.006211687295110 S10 0.014071338668160 S25 0.006449970233290 S9 0.015367550448190 S24 0.006710005953340 S8 0.016926489910360 S23 0.006991332142660 S7 0.018836924240150 S22 0.007297385745380 S6 0.021232615151970 S21 0.007631431941830 S5 0.024324905541030 S20 0.007997523135440 S4 0.028468352225130 S19 0.008400495372910 S4 0.034306329554970 S18 0.008846216714890 S3 0.043138734089010 S17 0.009341867768130 S2 0.058038919848870 S16 0.009896332328730 S1 0.088392216030230 S15 0.010520733534250

数据从I

29～I

1

单调减少，且S

1

和I

1

吻合得相当好，算法稳定，且精度高。

3.1由四条边四个初始点绘制图形

MATLAB程序：

clear

p=[0 0;0 10;10 10;10 0;0 0]; % the initial 4 points

n=4; % the initial 4 sides

A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)];

% matrix to cause 60 degree rotation

for k=1:4

j=0;

for i=1:n % i represents the number of sides during each cycling