《数值分析》课程设计—作业.doc

数值分析课程设计作业（本文档内的有些运行结果，限于篇幅，使文档结构更和谐、紧凑，已做相关的改动，程序代码没变） 实验一1.1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。解：一、 问题分析：对于本题，比较简单，我们只需要判断原来椰子的个数及每个人私藏了一份之后剩下的是否能被5除余1，直到最后分完。function fentao(n)a=cat(1,7);for j=n:-1:1 a(1)=j;i=1; while i<7 a(i+1)=4*(a(i)-1)/5; i=i+1; end if a(7)=fix(a(7) a, end endend二、 问题求解：通过matlab建立M文件，有如下程序：n=input('input n:');for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end if p=fix(p) break; endenddisp(x,p) 或者 input n:2000 1023 15621>> fentao(20001)a = 15621 12496 9996 7996 6396 5116 4092对于第一个程序，n取2000；对于第二个程序，n取20001，就能得到我们想要的结果，即原先一共有15621个椰子，最终平均每人得4092个椰子。 1.2 当时，选择稳定的算法计算积分.解：一、问题分析：由知： 以及： 得递推关系： ， 但是通过仔细观察就能知道上述递推公式每一步都将误差放大十倍，即使初始误差很小，但是误差的传播会逐步扩大，也就是说用它构造的算法是不稳定的，因此我们改进上述递推公式（算法）如下： 通过比较不难得出该误差是逐步缩小的，即算法是稳定的。二、 问题求解：为了利用上面稳定的算法，需要我们估计初值的值。因为 所以当n=100的时，我们有： 9.090909090909091e-0049.645173882220725e-004于是可取他们的平均值，有=9.368041486564908e-004，利用上述稳定算法，可求得相应的值和程序如下（限于篇幅，这里只给出了前后各十个连续的数据，详细的数据会连同作业一并上交）： 积分计算对照表 1n计算值准确值误差10.0591827610.059182761020.0402929470.040292947030.0294028930.0294028931.17997E-1640.0227087130.0227087131.5278E-1650.0183339560.018333956060.015312850.015312851.13285E-1670.0131250370.013125037080.01147650.0114765090.0101930630.0101930630910.000999990.000999992.68937E-11920.000989110.000989112.71896E-10930.0009784640.0009784642.74854E-09940.0009680450.0009680452.77813E-08950.0009578450.0009578462.80771E-07960.0009478620.0009478592.83729E-06970.0009380510.0009380782.86687E-05980.0009287660.0009284970.000289646990.0009164210.000919110.0029260391000.0009368040.0009099110.029556214 I=cat(1,100);J=cat(1,100);K=cat(1,100);I(100)=9.368041486564908e-004;format long;%求近似值for n=99:-1:1 I(n)=(1-exp(-n)/n-I(n+1)/10;end%求精确值for n=1:100 syms x; k=n; J(k)=int(exp(-k*x)/(exp(-x)+10),x,0,1);end%求误差for n=1:100 K(n)=abs(I(n)-J(n)/J(n);endn=1:100;A=n;I;J;K;B=A'xlswrite('1_2.xls',B)>> format long>> min=(1-exp(-100)/(11*100),max=(1-exp(-100)/(100*(exp(-1)+10)min = 9.090909090909091e-004max = 9.645173882220725e-004>> I100=(min+max)/2I100 = 9.368041486564908e-0041.3 绘制静态和动态的Koch分形曲线问题描述：从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形，这时，图形中共有17个结点。这种迭代继续进行下去可以形成Koch分形曲线。在迭代过程中，图形中的结点将越来越多，而曲线最终显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。Koch分形曲线的绘制与算法设计和计算机实现相关。图1.1 Koch曲线的形成过程解：一、算法分析： 考虑由直线段（2个点）产生第一个图形（5个点）的过程。上图中，设和分别为原始直线段的两个端点，现需要在直线段的中间依次插入三个点，。显然位于线段三分之一处，位于线段三分之二处，点的位置可看成是由点以点为轴心，逆时针旋转600而得。旋转由正交矩阵实现。算法根据初始数据（和点的坐标），产生图1中5个结点的坐标。结点的坐标数组形成一个矩阵，矩阵的第一行为的坐标，第二行为的坐标，第五行为的坐标。矩阵的第一列元素分别为5个结点的坐标，第二列元素分别为5个结点的坐标。进一步考虑该曲线形成过程中结点数目的变化规律。设第次迭代产生的结点数为，第次迭代产生的结点数为，则和中间的递推关系为。二、 问题求解: 、程序代码p=0 0;10 0; %P为初始两个点的坐标,第一列为x坐标,第二列为y坐标n=2; %n为结点数A=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); %旋转矩阵for k=1:6 d=diff(p)/3; %diff计算相邻两个点的坐标之差,得到相邻两点确定的向量 %则d就计算出每个向量长度的三分之一,与题中将线段三等分对应 m=4*n-3; %迭代公式 q=p(1:n-1,:); %以原点为起点,前n-1个点的坐标为终点形成向量 p(5:4:m,:)=p(2:n,:); %迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标 p(2:4:m,:)=q+d; %用向量方法计算迭代后处于4k+2位置上的点的坐标 p(3:4:m,:)=q+d+d*A' %用向量方法计算迭代后处于4k+3位置上的点的坐标 p(4:4:m,:)=q+2*d; %用向量方法计算迭代后处于4k位置上的点的坐标 n=m; subplot(2,3,k) plot(p(:,1),p(:,2),'r') %绘出每相邻两个点的连线axis(0,3,0,3) %迭代后新的结点数目 %axis offend 、运行结果： 实验二2.1 小行星轨道问题：一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在五个不同的对小行星作了五次观察，测得轨道上五个点的坐标数据（单位：万公里）如下表所示：P1P2P3P4P5X坐标5360558460628596666268894Y坐标602611179169542349268894由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：现需要建立椭圆的方程以供研究。(1) 分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组以及方程组的系数矩阵和右端项b；(2) 用MARLAB求低阶方程的指令Ab求出待定系数;(3) 分别用直接法、Jacobi迭代法、Gauss-Seidel迭代法求出待定系数.解：（1）、程序代码：X=53605 58460 62859 66662 68894;Y=6026 11197 16954 23492 68894;B=X.*X;2*X.*Y;Y.*Y;2*X;2*Y; %B是系数矩阵的转置A=B'%A是系数矩阵A=vpa(A,16)syms a1 a2 a3 a4 a5;a=a1;a2;a3;a4;a5;b=-1;-1;-1;-1;-1;b=sym(A*a) 、运行结果：A = 2873496025.0, 646047460.0, 36312676.0, 107210.0, 12052.0 3417571600.0, 1309153240.0, 125372809.0, 116920.0, 22394.0 3951253881.0, 2131422972.0, 287438116.0, 125718.0, 33908.0 4443822244.0, 3132047408.0, 551874064.0, 133324.0, 46984.0 4746383236.0, 9492766472.0, 4746383236.0, 137788.0, 137788.0 b = 2873496025.0*a1 + 646047460.0*a2 + 36312676.0*a3 + 107210.0*a4 + 12052.0*a5 3417571600.0*a1 + 1309153240.0*a2 + 125372809.0*a3 + 116920.0*a4 + 22394.0*a5 3951253881.0*a1 + 2131422972.0*a2 + 287438116.0*a3 + 125718.0*a4 + 33908.0*a5 4443822244.0*a1 + 3132047408.0*a2 + 551874064.0*a3 + 133324.0*a4 + 46984.0*a5 4746383236.0*a1 + 9492766472.0*a2 + 4746383236.0*a3 + 137788.0*a4 + 137788.0*a5 即是方程组的形式如下:-1=2873496025.0*a1 + 646047460.0*a2 + 36312676.0*a3 + 107210.0*a4 + 12052.0*a5-1=3417571600.0*a1 + 1309153240.0*a2 + 125372809.0*a3 + 116920.0*a4 + 22394.0*a5-1=3951253881.0*a1 + 2131422972.0*a2 + 287438116.0*a3 + 125718.0*a4 + 33908.0*a5-1=4443822244.0*a1 + 3132047408.0*a2 + 551874064.0*a3 + 133324.0*a4 + 46984.0*a5-1=4746383236.0*a1 + 9492766472.0*a2 + 4746383236.0*a3 + 137788.0*a4 + 137788.0*a5（2） 、用指令Ab求出待定系数的值如下：>> b=-1;-1;-1;-1;-1;A;a=Ab a = 0.000000000080151337137572070705034479362305 -0.000000000099919136582879577308894350837601 -0.00000000017639794583181340912897589334276 -0.000012556065177690492529375516352892 0.000015497775048603393791640753306649 （3） 、1）、直接法 、程序代码：function RA,RB,n,X=bilizy(A,b)B=A,b;n=length(b);RA=rank(A);RB=rank(B);rank_cha=RB-RA;if rank_cha>0, disp('因为RA!=RB，所以此方程组无解.') returnendif RA=RB if RA=n disp('因为RA=RB=n,所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); T=0;D=B's=max(abs(D(1:n,:); S=s' for p=1:n-1 Y,j=max(abs(B(p:n,p)./S(p:n); J=min(j); C=B(p,:);B(p,:)=B(J+p-1,:);B(J+p-1,:)=C; T=S(p);S(p)=S(J+p-1);S(J+p-1)=T; for k=p+1:n m=B(k,p)/B(p,p); B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1); A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); end else disp('请注意：因为RA=RB<n,所以此方程组有无穷多解.') endend 、运行结果：>>b=-1;-1;-1;-1;-1;RA,RB,n,X=bilizy(A,b)因为RA=RB=n,所以此方程组有唯一解.RA = 5RB = 5n = 5X = 1.0e-004 * 0.0000 -0.0000 -0.0000 -0.1256 0.1550>> X=vpa(X,16) X = 0.00000000008015133713756955 -0.0000000000999191365828799 -0.0000000001763979458318154 -0.00001255606517769044 0.00001549777504860351 2）、Jacobi迭代法 、程序代码：%P:范数的名称，P=1，2，inf%error:近似解a的误差%maxl:迭代的最大次数function a=Jacobi(A,b,a0,P,error,maxl)n n=size(A);a=zeros(n,1);for k=1:maxl for j=1:n a(j)=(b(j)-A(j,1:j-1,j+1:n) a0(1:j-1,j+1:n)/A(j,j); end erra=norm(a-a0,P); a0=a; if(erra<error) disp('近似解a分别是：') a=vpa(a,16); return endendif(erra>=error) disp('请注意：Jacobi迭代次数已经超过最大迭代次数maxl。')end 、运行结果：>>X=53605 58460 62859 66662 68894;Y=6026 11197 16954 23492 68894;B=X.*X;2*X.*Y;Y.*Y;2*X;2*Y; %B是系数矩阵的转置A=B'%A是系数矩阵b=-1;-1;-1;-1;-1;a0=0;0;0;0;0;>> a=Jacobi(A,b,a0,1,0.0001,100)近似解a分别是：a = -0.0000000003480081375786834 -0.0000000007638525188999265 -0.000000003479009721870011 -0.000007500525036752573 -0.000007257526054518536 3）、Gauss-Seidel迭代法 、程序代码：%P:范数的名称，P=1,2,inf;error：近似解x的误差;maxl：迭代的最大次数function a=GS(A,b,a0,P,error,maxl)n n=size(A);a=zeros(n,1);for k=1:maxl for j=1:n aa=0; for i=1:n if i<j aa=aa+A(j,i)*a(i); end if i>j aa=aa+A(j,i)*a0(i); end end a(j)=(b(i)-aa)/A(j,j); end erra=norm(a-a0,P); a0=a; if (erra<error) disp('近似解a分别是：') a=vpa(a,16); return endendif(erra>=error) disp('请注意：Gauss-Seidel迭代次数已经超过最大迭代次数maxl.')end 、运行结果：>> X=53605 58460 62859 66662 68894;Y=6026 11197 16954 23492 68894;B=X.*X;2*X.*Y;Y.*Y;2*X;2*Y; A=B'b=-1;-1;-1;-1;-1;a0=0;0;0;0;0;>> a=GS(A,b,a0,inf,0.00001,100)近似解a分别是：a = 0.0000000001678142444990526 0.000000002024408102712696 0.000000000369101766105352 -0.000009750490974607461 -0.00015547175272548952.3 设 (1) 将矩阵A进行LU分解A=LU，其中U是上三角矩阵，L是主对角线上的元素都是1的下三角矩阵。解：>> A=20 2 3;1 8 1;2 -3 15;>> L,U=lu(A)L = 1.0000 0 0 0.0500 1.0000 0 0.1000 -0.4051 1.0000U = 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443(2) 利用上述分解分别求解方程组并由此求出逆矩阵。>> b1=1;0;0;b2=0;1;0;b3=0;0;1;>> X1=U(Lb1), X2=U(Lb2),X3=U(Lb3),inv(A)X1 = 0.0517 -0.0055 -0.0080解：inv(A) = 0.0517 -0.0164 -0.0093 -0.0055 0.1237 -0.0072 -0.0080 0.0269 0.0665X3 = -0.0093 -0.0072 0.0665X2 = -0.0164 0.1237 0.0269(3) 用LU分解求下列线性方程组的解（方程组的精确解是） 解：>> A=4 2 -3 -1 2 1 0 0 0 0; 8 6 -5 -3 6 5 0 1 0 0; 4 2 -2 -1 3 2 -1 0 3 1; 0 -2 1 5 -1 3 -1 1 9 4; -4 2 6 -1 6 7 -3 3 2 3; 8 6 -8 5 7 17 2 6 -3 5; 0 2 -1 3 -4 2 5 3 0 1; 16 10 -11 -9 17 34 2 -1 2 2; 4 6 2 -7 13 9 2 0 12 4; 0 0 -1 8 -3 -24 -8 6 3 -1;>> b=5;12;3;2;3;46;13;38;19;-21;>> L U=lu(A);>> x=U(Lb)x = 1.0000 -1.0000 0.0000 1.0000 2.0000 0.0000 3.0000 1.0000 -1.0000 2.0000结果是：2.4 (1) 用追赶法求解方程组(a) ，(b) ，解： (a)、程序代码：%追赶法function x=zhuiganfa(A,b)n=rank(A);for i=1:n if A(i,i)=0 return; endendd=ones(n,1);a=ones(n-1,1);c=ones(n-1);for i=1:n-1 a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1)=A(n,n);%解Ly=b的解，解保存在b中for i=2:n d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1)*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1)*b(i-1,1);end%求解Uy=x的解x(n,1)=b(n,1)/d(n,1);for i=n-1:-1:1 x(i,1)=(b(i,1)-c(i,1)*x(i+1,1)/d(i,1);end 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1666 0.1668 0.1662 0.1683 0.1607 0.1891 0.0829 0.4793x = 0.4793 0.0829 0.1891 0.1607 0.1683 0.1662 0.1668 0.1666 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 运行结果：>> clear all>> A=zeros(30,30);A(30,30)=4;b=ones(30,1);b(1,1)=2;b(30,1)=2;for i=1:29 A(i,i+1)=1;A(i+1,i)=1; A(i,i)=4;end>> x=zhuiganfa(A,b) >> A=zeros(20,20);A(19,19)=5;A(20,19)=-2;A(19,20)=-2;A(20,20)=5;b=zeros(20,1);b(1,1)=1;for i=1:18 A(i,i+2)=1;A(i+2,i)=1; A(i,i+1)=-2;A(i+1,i)=-2; A(i,i)=5;end>> L,U=lu(A);>> x=vpa(U(Lb),6) 0.0000504376 0.000106309 0.0000292324 -0.0000143431 -0.0000126677 -0.00000146436 0.00000249126 0.00000131198 -9.33542*10(-8) -2.99738*10(-7)x = 0.242121 0.0943914 -0.0218242 -0.0313623 -0.00694256 0.00488693 0.00358612 0.000214823 -0.000784527 -0.000357862 (b)、 (2) 设计实验验证Hilbert矩阵的病态性，其中format rat;for n=1:10 A=hilb(n); tiaojianshuA=norm(A,1)*norm(inv(A),1)end解：tiaojianshuA = 1 tiaojianshuA = 27 tiaojianshuA = 748 tiaojianshuA = 28375 tiaojianshuA = 943656 tiaojianshuA = 29070279 tiaojianshuA =985194890 tiaojianshuA =33872791987 tiaojianshuA =1099651961846 tiaojianshuA =35353690673203 从上面的结果可以看出，Hilbert矩阵的条件数的增长速度惊人，才十阶，矩阵的条件数就为35353690673203，这是一个，可怕的数字，由它组成的方程组肯定是病态的。实验三3.1 用Taylor级数的第项多项式，分别在区间和上逼近正弦函数，并用计算机绘出上面31个函数的图形。解： 一、程序代码：syms x;y=sin(x);ezplot(y,0,pi);grid on;axis(0,pi,0,1.5);hold on;for m=1:2:30 p=taylor(y,x,m); ezplot(p,0,pi); axis(0,pi,0,1.5);end 二、函数图形： 3.2 已知四个点A，B，C，D的具体位置A(0, 0)，B(0, 3), C(8, 1), D(10, 5),求两个点H1，H2的具体位置，使AH1+BH1+H1H2+H2C+H2D为最短。解：function d =ju_li(a,b)d=sqrt(sum(a-b).2);returnendsyms x1 y1 x2 y2;A=0 0;B=0 3;C=8 1;D=10 5;H1=x1 y1;H2=x2 y2;sum_d=ju_li(A,H1)+ju_li(B,H1)+ju_li(H1,H2)+ju_li(C,H2)+ju_li(D,H2)x1 y1 x2 y2=solve(diff(sum_d,'x1'),diff(sum_d,'y1'),diff(sum_d,'x2'),diff(sum_d,'y2'),x1,y1,x2,y2)sum_d =3*x12 - 2*x1*x2 + 3*x22 - 36*x2 + 3*y12 - 2*y1*y2 - 6*y1 + 3*y22 - 12*y2 + 199x1 =9/4y1 =27/4x2 =15/8y2 =21/8实验四4.1 方程的根实际上是两个函数的交点的横坐标，用计算机绘出两个函数在区间的图形，观察图形，分析它们的交点分布规律及特点，试写出方程的全部实根所在的隔根区间；并用二分法求出每一个根的近似值。 解： 一、画函数图像如下：>> syms x y1 y2;>> y1=sin(x);y2=1/x;>> ezplot(y1,-6,6);grid on; hold on; ezplot(y2,-6,6)>> gtext('y1=sin(x)')>> gtext('y2=1/x') 二、 用二分法求近似解: 通过上面的函数图像可以知道，y1=sin(x)与y2=1/x在区间-6,6