西安交通大学2015年计算方法上机作业(共30页).docx

《西安交通大学2015年计算方法上机作业(共30页).docx》由会员分享，可在线阅读，更多相关《西安交通大学2015年计算方法上机作业(共30页).docx（30页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上计算方法上机作业学院：能动学院班级：硕4011班P81页第四题求20阶三对角方程组的解，其中：，算法分析：追赶法的算法组织如下：(1) 输入三对角矩阵和右端向量；(2) 将压缩为四个一维数组，将分解矩阵压缩为三个一维数组(3) 对做分解（也可以用分解）导出追赶法的计算步骤如下：(4) 回代求解 (5) 停止，输出结果 结果：将结果带入方程可以验证结果完全正确。第五题算法分析：将原始线性代数方程组改写为的形式，其中为的矩阵函数。于是可以得到迭代格式：，此即为迭代法的迭代格式。如果在计算时，将已经算出的分量立即代换对应分量，则得到迭代法的迭代格式。1. 迭代法的算法组织如

2、下：(1) 给出迭代格式(2) 给出迭代初始向量、允许误差和最大迭代次数(3) 按照迭代格式进行迭代，直至达满足迭代停止条件(4) 停止，输出结果2. 迭代法的算法组织如下：(1) 给出迭代格式(2) 给出迭代初始向量、允许误差和最大迭代次数(3) 按照迭代格式，并且将已经算出的分量立即代换对应分量进行迭代，直至达满足迭代停止条件(4) 停止，输出结果结果：将结果带入方程验证的确满足误差1e-3P131页第一题已知，对a 计算函数在点处的值b 求插值数据点的插值多项式和三次样条插值多项式c 对，计算和相应的函数值d 计算，解释所得到结果。算法分析1. 插值插值多项式是实质是插值多项式的零次式和

3、一次式的推广，一般形式为：2. 三次样条插值在节点分成的每个小区间上利用其节点处的二次导数值（即弯矩值）进行线性插值，再在此区间上积分两次，利用节点处的函数值求得两个积分常数，从而得到三次样条函数的表达式。要确定三次样条插值函数，用到的光滑性，利用三次样条函数在节点处的一阶导数连续得到满足个方程：其中三弯矩方程组只有个方程，不能确定未知量。为此，可以给出区间的端点处的一阶导数值。由和可以导出其中联合先前得到的三弯矩方程组和两个补充方程得到恰定的三弯矩方程组。它是严格对角占优的三对角方程组，可以用追赶法求解。结果：a 直接计算即可得出对于n=5有对于n=10有对于n=20有b对于n=5有牛顿插值

4、多项式样条插值函数（分别对于各个小区间）为：对于n=10有牛顿插值函数p = (84*X9 + *X8 - 112*X7 - *X6 + 136*X5 + *X4 + 44*X3 - *X2 + 4*X + 93257)/26848 pp = 3.73*10(-14)*X9 + 21.6*X8 - 6.95*10(-14)*X7 - 44.9*X6 + 4.16*10(-14)*X5 + 30.7*X4 + 7.86*10(-15)*X3 - 8.26*X2 + 1.32*10(-15)*X + 0.862对于n=10有样条插值函数：对于n=20有牛顿插值函数- 6.82*10(-10)*X1

5、9 - 25770.0*X18 + 2.52*10(-9)*X17 + 95660.0*X16 - 3.79*10(-9)*X15 - 1.47*105*X14 + 3.02*10(-9)*X13 + 1.21*105*X12 - 1.46*10(-9)*X11 - 58660.0*X10 + 4.3*10(-10)*X9 + 17220.0*X8 - 6.6*10(-11)*X7 - 3066.0*X6 + 4.06*10(-12)*X5 + 328.0*X4 - 3.19*10(-13)*X3 - 21.6*X2 + 3.08*10(-15)*X + 0.993对于n=20有样条插值函数：

6、 - 5.25*X3 - 14.37*X2 - 12.98*X - 3.827 1.723*X3 + 4.461*X2 + 3.963*X + 1.256 0.07504*X3 + 0.5066*X2 + 0.7996*X + 0.4127 0.9738*X3 + 2.394*X2 + 2.121*X + 0.721 1.552*X3 + 3.436*X2 + 2.746*X + 0.846 3.552*X3 + 6.434*X2 + 4.245*X + 1.096 5.726*X3 + 9.044*X2 + 5.289*X + 1.235 12.54*X3 + 15.17*X2 + 7.12

7、7*X + 1.419 - 32.79*X3 - 12.02*X2 + 1.688*X + 1.056 - 89.07*X3 - 28.91*X2 - 2.345*10(-15)*X + 1.0 89.07*X3 - 28.91*X2 - 6.395*10(-15)*X + 1.0 32.79*X3 - 12.02*X2 - 1.688*X + 1.056 - 12.54*X3 + 15.17*X2 - 7.127*X + 1.419 - 5.726*X3 + 9.044*X2 - 5.289*X + 1.235 - 3.552*X3 + 6.434*X2 - 4.245*X + 1.096

8、- 1.552*X3 + 3.436*X2 - 2.746*X + 0.846 - 0.9738*X3 + 2.394*X2 - 2.121*X + 0.721 - 0.07504*X3 + 0.5066*X2 - 0.7996*X + 0.4127 - 1.723*X3 + 4.461*X2 - 3.963*X + 1.256 5.25*X3 - 14.37*X2 + 12.98*X - 3.827c 容易计算各点函数值 YK = 0.0400 0.0416 0.0433 0.0451 0.0471 0.0491 0.0513 0.0536 0.0561 0.0588 0.0617 0.06

9、48 0.0681 0.0716 0.0755 0.0796 0.0841 0.0890 0.0943 0.1000 0.1063 0.1131 0.1206 0.1289 0.1379 0.1479 0.1590 0.1712 0.1848 0.2000 0.2169 0.2358 0.2571 0.2809 0.3077 0.3378 0.3717 0.4098 0.4525 0.5000 0.5525 0.6098 0.6711 0.7353 0.8000 0.8621 0.9174 0.9615 0.9901 1.0000 0.9901 0.9615 0.9174 0.8621 0.8

10、000 0.7353 0.6711 0.6098 0.5525 0.5000 0.4525 0.4098 0.3717 0.3378 0.3077 0.2809 0.2571 0.2358 0.2169 0.2000 0.1848 0.1712 0.1590 0.1479 0.1379 0.1289 0.1206 0.1131 0.1063 0.1000 0.0943 0.0890 0.0841 0.0796 0.0755 0.0716 0.0681 0.0648 0.0617 0.0588 0.0561 0.0536 0.0513 0.0491 0.0471 0.0451 0.0433 0.

11、0416 0.0400 对于n=5 可以计算相应插值函数的值Nk= -0.0496 -0.1257 -0.1906 -0.2449 -0.2891 -0.3239 -0.3497 -0.3672 -0.3769 -0.3793 -0.3750 -0.3643 -0.3479 -0.3263 -0.2997 -0.2688 -0.2340 -0.1957 -0.1542 -0.1101 -0.0636 -0.0152 0.0347 0.0859 0.1379 0.1905 0.2434 0.2962 0.3487 0.4005 0.4515 0.5014 0.5499 0.5968 0.6419

12、 0.6850 0.7260 0.7646 0.8008 0.8342 0.8649 0.8927 0.9174 0.9391 0.9576 0.9728 0.9846 0.9932 0.9983 1.0000 0.9983 0.9932 0.9846 0.9728 0.9576 0.9391 0.9174 0.8927 0.8649 0.8342 0.8008 0.7646 0.7260 0.6850 0.6419 0.5968 0.5499 0.5014 0.4515 0.4005 0.3487 0.2962 0.2434 0.1905 0.1379 0.0859 0.0347 -0.01

13、52 -0.0636 -0.1101 -0.1542 -0.1957 -0.2340 -0.2688 -0.2997 -0.3263 -0.3479 -0.3643 -0.3750 -0.3793 -0.3769 -0.3672 -0.3497 -0.3239 -0.2891 -0.2449 -0.1906 -0.1257 -0.0496 Sk=0.0383 0.0378 0.0371 0.0363 0.0355 0.0348 0.0344 0.0343 0.0347 0.0357 0.0373 0.0397 0.0430 0.0472 0.0526 0.0592 0.0672 0.0765

14、0.0874 0.1000 0.1143 0.1302 0.1475 0.1660 0.1857 0.2063 0.2277 0.2497 0.2721 0.2948 0.3175 0.3403 0.3628 0.3849 0.4064 0.4273 0.4472 0.4661 0.4837 0.5000 0.5147 0.5279 0.5396 0.5497 0.5582 0.5652 0.5706 0.5745 0.5768 0.5776 0.5768 0.5745 0.5706 0.5652 0.5582 0.5497 0.5396 0.5279 0.5147 0.5000 0.4837

15、 0.4661 0.4472 0.4273 0.4064 0.3849 0.3628 0.3403 0.3175 0.2948 0.2721 0.2497 0.2277 0.2063 0.1857 0.1660 0.1475 0.1302 0.1143 0.1000 0.0874 0.0765 0.0672 0.0592 0.0526 0.0472 0.0430 0.0397 0.0373 0.0357 0.0347 0.0343 0.0344 0.0348 0.0355 0.0363 0.0371 0.0378 0.0383 对于n=20 可以计算相应插值函数的值Nk= 8.2065 7.7

16、749 4.8826 2.1401 0.3516 -0.4925 -0.6818 -0.5282 -0.2658 -0.0308 0.1204 0.1843 0.1839 0.1492 0.1065 0.0729 0.0560 0.0558 0.0680 0.0865 0.1057 0.1218 0.1334 0.1407 0.1454 0.1497 0.1557 0.1648 0.1779 0.1948 0.2149 0.2373 0.2611 0.2861 0.3123 0.3403 0.3713 0.4066 0.4476 0.4950 0.5492 0.6094 0.6741 0.74

17、08 0.8063 0.8669 0.9189 0.9589 0.9841 0.9927 0.9841 0.9589 0.9189 0.8669 0.8063 0.7408 0.6741 0.6094 0.5492 0.4950 0.4476 0.4066 0.3713 0.3403 0.3123 0.2861 0.2611 0.2373 0.2149 0.1948 0.1779 0.1648 0.1557 0.1497 0.1454 0.1407 0.1334 0.1218 0.1057 0.0865 0.0680 0.0558 0.0560 0.0729 0.1065 0.1492 0.1

18、839 0.1843 0.1204 -0.0308 -0.2658 -0.5282 -0.6818 -0.4925 0.3516 2.1401 4.8826 7.7749 8.2065 Sk=0.0390 0.0403 0.0423 0.0446 0.0471 0.0494 0.0516 0.0539 0.0563 0.0588 0.0616 0.0647 0.0680 0.0716 0.0755 0.0796 0.0841 0.0890 0.0943 0.1000 0.1063 0.1131 0.1206 0.1289 0.1379 0.1479 0.1589 0.1712 0.1848 0

19、.2000 0.2170 0.2359 0.2572 0.2810 0.3077 0.3375 0.3711 0.4089 0.4517 0.5000 0.5541 0.6127 0.6742 0.7372 0.8000 0.8606 0.9152 0.9594 0.9891 1.0000 0.9891 0.9594 0.9152 0.8606 0.8000 0.7372 0.6742 0.6127 0.5541 0.5000 0.4517 0.4089 0.3711 0.3375 0.3077 0.2810 0.2572 0.2359 0.2170 0.2000 0.1848 0.1712

20、0.1589 0.1479 0.1379 0.1289 0.1206 0.1131 0.1063 0.1000 0.0943 0.0890 0.0841 0.0796 0.0755 0.0716 0.0680 0.0647 0.0616 0.0588 0.0563 0.0539 0.0516 0.0494 0.0471 0.0446 0.0423 0.0403 0.0390 d 有c得出的结果很容易计算对于n=5个节点时E（N）= 0.4381 E (s)= 0.4224对于n=20个节点时E（N）= 8.1665E (s)= 0.0031这是由于已知数据的误差对插值多项式计算的结果带来了很大

21、的误差。具体说，当函数值有误差而变为且误差限为。由于插值与插值本质上是一样的，则插值与插值的误差估计是一样的。由插值多项式的计算值也会有误差，于是可以导出插值多项式的误差当很大时，的误差就很大，于是就出现了前面观察到的现象。p189第三题 算法分析：此方法称为复化积分法那么有了算法就可以计算了。结果：利用simpson公式可以计算得出如下的方程式，求解方程式便可得出相应函数节点上的值。y1为函数y（x）在x=0处的值y2为函数y（x）在x=0.125处的值以此类推y9为函数y（x）在x=1处的值可以得出下列方程(1) y1/24 + (4*y2)/27 + y3/15 + (4*y4)/33

22、+ y5/18 + (4*y6)/39 + y7/21 + (4*y8)/45 + y9/48 + 5/8=y1 (2) (7*y1)/192 + (55*y2)/432 + (9*y3)/160 + (53*y4)/528 + (13*y5)/288 + (17*y6)/208 + (25*y7)/672 + (49*y8)/720 + y9/64 + 619/1296=y2 (3) y1/32 + (23*y2)/216 + (11*y3)/240 + (7*y4)/88 + (5*y5)/144 + (19*y6)/312 + (3*y7)/112 + (17*y8)/360 + y9/

23、96 + 39/100=y3 (4)(5*y1)/192 + (37*y2)/432 + (17*y3)/480 + (31*y4)/528 + (7*y5)/288 + (25*y6)/624 + (11*y7)/672 + (19*y8)/720 + y9/192 + 661/1936=y4 (5) y1/48 + (7*y2)/108 + y3/40 + (5*y4)/132 + y5/72 + y6/52 + y7/168 + y8/180 + 23/72=y5 (6) y1/64 + (19*y2)/432 + (7*y3)/480 + (3*y4)/176 + y5/288 - y

24、6/624 - y7/224 - (11*y8)/720 - y9/192 + 855/2704=y6 (7) y1/96 + (5*y2)/216 + y3/240 - y4/264 - y5/144 - (7*y6)/312 - (5*y7)/336 - (13*y8)/360 - y9/96 + 16/49=y7 (8) y1/192 + y2/432 - y3/160 - (13*y4)/528 - (5*y5)/288 - (9*y6)/208 - (17*y7)/672 - (41*y8)/720 - y9/64 + 1249/3600=y8 (9) 3/8 - y3/60 - y

25、4/22 - y5/36 - (5*y6)/78 - y7/28 - (7*y8)/90 - y9/48 - y2/54=y9容易解得函数f（x）在各个节点处的值用三次样条函数近似真解90.*X3 - 24.*X2 - 0.15625*X + 1.26.*X3 + 19.*X2 - 5.*X + 1.5.*X3 - 5.*X2 + 0.*X + 0.2.*X3 + 4.*X2 - 2.*X + 1.0.*X3 + 0.*X2 - 1.*X + 0. 0.*X3 - 0.*X2 - 0.*X + 0. - 3.*X3 + 7.*X2 - 6.*X + 2. 11.*X3 - 30.*X2 + 2

26、7.*X - 7.可以计算插值函数在节点处的值为1.0000 0.5729 0.0972 0.8056 0.9897 0.0000 0.0309 0.5552 0.1811 0.2828 0.0000 0.86010.6760 0.7305 0.96910.0000 0.0103 0.90870.0457 0.2842 0.0000真解在节点处的值为 1.0000 0.8050 0.1735 0.6787 0.4444 0.0000 0.1361 0.3141 0.2653 0.7099 0.4444 0.0697 0.0000 0.8549 0.5675 0.4898 0.8642 0.78

27、60 0.9308 0.9494 0.0000绿色为真解红色为插值解可以看出他们非常相近p229 第三题用Newton迭代法求解方程组算法分析：用迭代法解非线性方程时，如何构造迭代函数是非常重要的，那么怎样构造的迭代函数才能保证迭代法收敛呢？牛顿迭代法就是常用的方法之一，其迭代格式的来源大概有以下几种方式：1设，对在点作泰勒展开：略去二次项，得到的线性近似式：。由此得到方程0的近似根(假定0)，即可构造出迭代格式(假定0)： 这就是牛顿迭代公式，若得到的序列收敛于，则就是非线性方程的根。2 牛顿迭代法也称为牛顿切线法，这是由于的线性化近似函数是曲线过点的切线而得名的，求的零点代之以求的零点，即

28、切线与轴交点的横坐标，如右图所示，这就是牛顿切线法的几何解释。实际上，牛顿迭代法也可以从几何意义上推出。利用牛顿迭代公式，由得到，从几何图形上看，就是过点作函数的切线，切线与轴的交点就是，所以有，整理后也能得出牛顿迭代公式：。3 要保证迭代法收敛，不管非线性方程0的形式如何，总可以构造： 作为方程求解的迭代函数。因为：而且在根附近越小，其局部收敛速度越快，故可令：若0(即根不是0的重根)，则由得：，因此可令，则也可以得出迭代公式：。4 迭代法的基本思想是将方程改写成等价的迭代形式，但随之而来的问题却是迭代公式不一定收敛，或者收敛的速度较慢。运用前述加速技巧，对于简单迭代过程，其加速公式具有形式

29、：，其中记，上面两式可以合并写成：这种迭代公式称作简单的牛顿公式，其相应的迭代函数是：。需要注意的是，由于是的估计值，若取，则实际上便是的估计值。假设，则可以用代替上式中的，就可得到牛顿法的迭代公式：。牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程逐步归结为某种线性方程来求解。对于方程组的解法和单个方程的解法类似。第一个方程计算结果：format long 格式后更精确地解如下将所算结果代入方程可以得到误差满足程序设定的误差要求。第二个方程的计算结果：format long 格式后更精确地解如下第287页第8题用Runge-Kutta方法 计算初值问题 取步长h=0.25算法分析：

30、对于一阶常微分方程初值问题 （1）的数值解法是近似计算中很中的一部分。常微分方程的数值解法通常就是给出定义域上的n个等距节点，求出所对应的函数值yn。通常其数值解法可分为两大类：1. 单步法：这类方法在计算yn+1的值时，只需要知道xn+1、xn和yn即可，就可算出。这类方法典型有欧拉法和R-K法。2. 多步法：这类方法在计算yn+1的值时，除了需要知道xn+1、xn和yn值外，还需要知道前k步的值。典型的方法如Adams法。经典的R-K法是一个四阶的方法。它最大的优点就是它是单步法，精度高，计算过程便于改变步长。其缺点也很明显，计算量大，每前进一步就要计算四次函数值f。它的具体的计算公式如下

31、式所示。 结果：不妨减小步长，算出更为精确地结果 h=0.01时 结果如下图function Gaussseidel% 81页第五题用Gauss-seidel迭代求解电路电流方程组，使误差小于1e-3clcclear allclose allA=28 -3 0 0 0 ;-3 38 -10 0 -5; 0 -10 25 -15 0; 0 0 -15 45 0;0 -5 0 0 30b=10 0 0 0 0ep=1e-3;x=maseidel(A,b,ep)function x=maseidel(A,b,ep,x0)% x0为初始迭代值，ep为最终解得误差n=length(b);if nargi

32、n4;%若输入变量小于4个默认初始迭代值为零向量。 x0=zeros(n,1);endx=zeros(n,1);k=0;while 1 for i=1:n if i=1 x(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1); else if i=n x(n)=(b(n)-A(n,1:n-1)*x(1:n-1)/A(n,n); else x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i); end end end if norm(x-x0,inf)ep break end x0=x; k=k+1;endfuncti

33、on Jacobi% 81页第五题用Jacobi迭代求解电路电流方程组，使误差小于1e-3clear allclose allclcA=28 -3 0 0 0 ;-3 38 -10 0 -5; 0 -10 25 -15 0; 0 0 -15 45 0;0 -5 0 0 30b=10 0 0 0 0ep=1e-3;x=majacobi(A,b,ep)function x=majacobi(A,b,ep,x0)% x0为初始迭代值，ep为最终解得误差n=length(b);if nargin4;% 若输入变量小于4个默认初始迭代值为零向量。 x0=zeros(n,1);endx=zeros(n,1

34、);k=0;while 1 for i=1:n x(i)=(b(i)-A(i,1:i-1,i+1:n)*x0(1:i-1,i+1:n)/A(i,i); end if norm(x-x0,inf)ep%判断迭代前与迭代后之间的误差是否小于要求的误差 break end x0=x; k=k+1;endfunction zhuiganfaclcclear allclose all% 81页 第四题解法a=-1*ones(20,1);b=4*ones(20,1);c=-1*ones(20,1);d=2*ones(20,1);d(1)=3;d(20)=3;x=machase(a,b,c,d)functi

35、on x=machase(a,b,c,d)% 用途:追赶法解三对角方程组Ax=d% 格式：x=machase(a,b,c,d) a为次下对角线元素向量，b%为主对角线元素向量，c为次上对角线元素向量，d为右端向量，x返回解向量n=length(a);%追for k=2:n b(k)=b(k)-a(k)*c(k-1)/b(k-1); d(k)=d(k)-a(k)*d(k-1)/b(k-1);end%赶x(n)=d(n)/b(n);for k=n-1:-1:1 x(k)=(d(k)-c(k)*x(k+1)/b(k);end clc;clear all;close all;n=4;h=1/n;%步长

36、为0.25for i=1:n+1 t(i)=(i-1)*h;%0t 1endf1 = (t, y)(1-t.2-y.2)y(1)=1; for i=1:n K1=h*f1(t(i),y(i); K2=h*f1(t(i)+0.5*h,y(i)+0.5*K1); K3=h*f1(t(i)+0.5*h,y(i)+0.5*K2); K4=h*f1(t(i)+h,y(i)+K3); y(i+1)=y(i)+1/6*(K1+2*K2+2*K3+K4);endplot(t,y,b) tytitle(Runge-Kutta 算法计算结果) xlabel(T axis)ylabel(Y axis)clearcl

37、cclose all%程序计算牛顿插值多项式f=(x)1/(1+25*x.2)a=-1;b=1;n=20;%n=10， n=20 分别修改n值然后计算for i=1:nx(i)=a+(b-a)/(n-1)*(i-1);y(i,1)=eval(subs(f,x,x(i);endyfor j=1:n-1for k=j:n-1temp=y(k+1,j)-y(k,j);y(k+1,j+1)=temp/(x(k+1)-x(k+1-j);endc(j)=y(j,j);c(j+1)=y(j+1,j+1);endp=c(1);q=1;syms Xfor i=2:n q=q*(X-x(i-1); p=p+c(i

38、)*q;endp=simple(p)pp=vpa(p,3)for i=1:101t(i)=a+(b-a)/100*(i-1);Nn(i)=eval(subs(p,X,t(i);endfor i=1:101h(i)=a+(b-a)/100*(i-1);yy(i)=eval(subs(f,x,h(i);endNn% 画出各个图像plot(h,yy,r)hold onplot(t,Nn,b)hold ongrid onxlabel(x)ylabel(f(x)function yi=spline(x,y,ydot,xi)% ydot为给定的左右两个区间样条插值函数的一阶导数值 % p为所求样条插值函数，xi为插值节点，yi为插值函数在插值节点的值format longn=length(x);h=zeros(1,n);lambda=ones(1,n);mu=ones(1,n);M=zeros(n,1);d=zeros(n,1);for k=2:n h(k)=x(k)-x(k-1);endfor k=2:n-1 lambda(k)=h(k+1)/(h(k)+h(k+1); mu(k)=1-lambda(k); d(k)=6/( h(k)+h(k+1) )*( (y(k+1)-y(k) /h(k+1)-(y(k)-y(k-1)/h(k);endd(1)=6/h(2)*(