数值分析上机实验.doc

刘力辉 2010210804011.“画圆为方”问题也是古希腊人所提出几何三大难题中的另一个问题。即求作一个正方形，使其面积等于已知圆的面积。不妨设已知圆的半径为 R = 1，试用数值试验显示“画圆为方”问题计算过程中的误差。（1）MATLAB程序：y=pi(1/2); % to generate 15-bit value of square root of pib=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)=x3; l(k)=x;endformat longl,s（2） 误差分析：位数 h近似值V近似值2 1.7 4.913000000000003 1.775.545233000000004 1.7725.564051648000005 1.77245.567820479424006 1.772455.568291702981137 1.7724535.568319977240018 1.77245385.568327517058549 1.772453855.568327988297422 算法的数值稳定性实验设， 由 xn = xn + 5 xn 1 5 xn 1 可得递推式In = 5In 1 + 1/ n （1）从 I0 尽可能精确的近似值出发，利用递推公式： ( n = 1,2,20)计算从 I1 到 I20 的近似值； （2）从 较粗糙的估计值出发，用递推公式： ( n= 30,29, , 3, 2 )计算从 到 的近似值；(3) 分析所得结果的可靠性以及出现这种现象的原因。I0 =ln(5+x)=ln6-ln5 所以I00.18232155679395MATLAB 程序：format long I0=log2(6)/log2(exp(1)-log2(5)/log2(exp(1) % calculate the value of I0=ln6-ln5for n=1:20 I0=-5*I0+1/n; % recycling equation between I(n+1) and I(n) s(n)=I0;ends则计算结果为：表1I10.0883922160302300 I110.0140713362538500 I20.0580389198488700 I120.0129766520640700 I30.0431387340890000 I130.0120398166027400 I40.0343063295550100 I140.0112294884148600 I50.0284683522249700 I150.0105192245923700 I60.0243249055418100 I160.0099038770381400 I70.0212326151481100 I170.0093041442210800 I80.0188369242594600 I180.0090348344501700 I90.0169264898137900 I190.0074574066965100 I100.0153675509310500 I200.0127129665174600 从计算的数据看出I20=0.0127129665174600 I19=0.0074574066965100又In的积分范围为01，所以应该有InIn+1。所以算法不稳定。下面分析导致算法不稳定的原因：令Sn为近似值，则有 (1)又 (2)(1) -(2)得Sn-In=-5Sn-1-In-1=(-5)n(S0-I0)=(-5)n0 (3)0为I0的误差。可知由下往上递推时误差是以指数增加传递的，即n=(-5)n0 所以必然会导致算法不稳定。又公式(3)知，当由上往下递推时，误差是以指数减少传递的，即0=1/(-5)nn (4)所以此时算法是稳定的。下面便给出由I30的估计值算I29I1的值。 所以有In 即1/186I301/155，取为1/155。MATLAB程序：format long I0=1/155 % calculate the value of I30=ln6-ln5for n=30:-1:2 I0=(-1/5)*I0+1/(5*n); % recycling equation between I(n+1) and I(n) s(n-1)=I0;endj=1;for i=29:-1:1 a(j)=s(i);j=j+1;enda计算出相应数据为：表2S290.005376344086020 S140.011229186626480 S280.005821282906930 S130.012039876960420 S270.005978600561470 S120.012976639992530 S260.006211687295110 S100.014071338668160 S250.006449970233290 S90.015367550448190 S240.006710005953340 S80.016926489910360 S230.006991332142660 S70.018836924240150 S220.007297385745380 S60.021232615151970 S210.007631431941830 S50.024324905541030 S200.007997523135440 S40.028468352225130 S190.008400495372910 S40.034306329554970 S180.008846216714890 S30.043138734089010 S170.009341867768130 S20.058038919848870 S160.009896332328730 S10.088392216030230 S150.010520733534250 数据从I29I1 单调减少，且S1和I1吻合得相当好，算法稳定，且精度高。3. 1由四条边四个初始点绘制图形MATLAB程序：clearp=0 0;0 10;10 10;10 0;0 0; % the initial 4 pointsn=4; % the initial 4 sidesA=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); % matrix to cause 60 degree rotationfor k=1:4 j=0; for i=1:n % i represents the number of sides during each cycling q1=p(i,:); q2=p(i+1,:); d=(q2-q1)/3; j=j+1;r(j,:)=q1; % j represents the number of points the sides generate during each cycling j=j+1;r(j,:)=q1+d; j=j+1;r(j,:)=q1+d+d*A; j=j+1;r(j,:)=q1+2*d; end n=4*n;clear p p=r;q2;endplot(p(:,1),p(:,2)生成图片：图13.2 MATLAB程序：clearp=0 0;10 0; % the initial 4 pointsn=1; % the initial 4 sidesA1=cos(pi/2) -sin(pi/2);sin(pi/2) cos(pi/2); % matrix to cause 60 degree ratationA2=cos(pi/4) -sin(pi/4);sin(pi/4) cos(pi/4);for k=1:5 j=0; for i=1:n % i represents the number of sides during each cycling q1=p(i,:); q2=p(i+1,:); d=(q2-q1)/3; j=j+1;r(j,:)=q1; % j represents the number of points the sides generate during each cycling j=j+1;r(j,:)=q1+d; j=j+1;r(j,:)=q1+d+d*A1; j=j+1;r(j,:)=q1+d+2(1/2)*d*A2; j=j+1;r(j,:)=q1+2*d; end n=5*n;clear p p=r;q2endplot(p(:,1),p(:,2)生成图片：图2感觉图2没有达到预期的效果,这是怎么回事呢？仔细看看原来是边取的等距离即d一样，在以后生成的正方形就会发生交叠。所以为了避免这种情况发生，线段不取三等分。修干如下：程序：clearp=0 0;10 0; % the initial 4 pointsn=1; % the initial 4 sidesA1=cos(pi/2) -sin(pi/2);sin(pi/2) cos(pi/2); % matrix to cause 60 degree ratationA2=cos(pi/4) -sin(pi/4);sin(pi/4) cos(pi/4);for k=1:5 j=0; for i=1:n % i represents the number of sides during each cycling q1=p(i,:); q2=p(i+1,:); d=(q2-q1)/5; j=j+1;r(j,:)=q1; % j represents the number of points the sides generate during each cycling j=j+1;r(j,:)=q1+2*d; j=j+1;r(j,:)=q1+2*d+d*A1; j=j+1;r(j,:)=q1+2*d+2(1/2)*d*A2; j=j+1;r(j,:)=q1+3*d; end n=5*n;clear p p=r;q2endplot(p(:,1),p(:,2)生成图片：图3生成的图片就跟图1样式差不多了。3.3以下图形为生成元MATLAB程序：clearp=0 0;0 10;10 10;10 0;0 0; % the initial 4 pointsn=4; % the initial 4 sidesA=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); % matrix to cause 60 degree rotation clolkwisefor k=1:6 j=0; for i=1:n % i represents the number of sides during each cycling q1=p(i,:); q2=p(i+1,:); d=(q2-q1)/3; j=j+1;r(j,:)=q1; % j represents the number of points the sides generate during each cycling j=j+1;r(j,:)=q1+d; j=j+1;r(j,:)=q1+2*d+d*A; j=j+1;r(j,:)=q1+2*d; end n=4*n;clear p p=r;q2;endplot(p(:,1),p(:,2)生成图片：图43.4生成元图形如下：MATLAB程序：p=0 0;10 0; %P为初始两个点的坐标n=2; %n为结点数A=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); B=cos(-pi/3) -sin(-pi/3);sin(-pi/3) cos(-pi/3); C=cos(-pi/2) -sin(-pi/2);sin(-pi/2) cos(-pi/2);%A使向量顺时针旋转60度，B使向量逆时针旋转60度，C使向量顺时针旋转90度for k=1:4 d=diff(p)/3; d1=d(1:2:n,:);%取每条线段对应的向量 m=5*n; %迭代公式，即新生成点数 q1=p(1:2:n-1,:); p(10:10:m,:)=p(2:2:n,:); p(1:10:m,:)=p(1:2:n,:); %迭代后处于10k与10k+1位置上的点的坐标为迭代前的相应坐标p(2:10:m,:)=q1+d1; %用向量方法计算迭代后处于10k+2,10k+3,10k+5位置上的点的坐标,都相同 p(3:10:m,:)=p(2:10:m,:); p(4:10:m,:)=q1+d1+d1*A; %用向量方法计算迭代后处于10k+4位置上的点的坐标 p(5:10:m,:)=p(2:10:m,:);p(6:10:m,:)=q1+2*d1; %用向量方法计算迭代后处于10k+6,10k+7,10k+9位置上的点的坐标,都