西安交通大学2015年计算方法上机作业.doc

计算方法上机作业

学院：能动学院
班级：硕4011班

P81页第四题
求20阶三对角方程组的解，其中：
，
算法分析：
追赶法的算法组织如下：
(1) 输入三对角矩阵和右端向量；
(2) 将压缩为四个一维数组，将分解矩阵压缩为三个一维数组
(3) 对做分解（也可以用分解）导出追赶法的计算步骤如下：

(4) 回代求解

(5) 停止，输出结果

结果：

将结果带入方程可以验证结果完全正确。

第五题

算法分析：
将原始线性代数方程组改写为的形式，其中为的矩阵函数。于是可以得到迭代格式：，此即为迭代法的迭代格式。如果在计算时，将已经算出的分量立即代换对应分量，则得到迭代法的迭代格式。
1. 迭代法的算法组织如下：
(1) 给出迭代格式
(2) 给出迭代初始向量、允许误差和最大迭代次数
(3) 按照迭代格式进行迭代，直至达满足迭代停止条件
(4) 停止，输出结果
2. 迭代法的算法组织如下：
(1) 给出迭代格式
(2) 给出迭代初始向量、允许误差和最大迭代次数
(3) 按照迭代格式，并且将已经算出的分量立即代换对应分量进行迭代，直至达满足迭代停止条件
(4) 停止，输出结果

结果：

将结果带入方程验证的确满足误差<1e-3

P131页第一题
已知，对

a 计算函数在点处的值
b 求插值数据点的插值多项式和三次样条插值多项式
c 对，计算和相应的函数值
d 计算，，解释所得到结果。
算法分析
1. 插值
插值多项式是实质是插值多项式的零次式和一次式的推广，一般形式为：

2. 三次样条插值
在节点分成的每个小区间上利用其节点处的二次导数值（即弯矩值）进行线性插值，再在此区间上积分两次，利用节点处的函数值求得两个积分常数，从而得到三次样条函数的表达式。
要确定三次样条插值函数，用到的光滑性，利用三次样条函数在节点处的一阶导数连续得到满足个方程：

其中
三弯矩方程组只有个方程，不能确定未知量。为此，可
以给出区间的端点处的一阶导数值。
由和可以导出

其中
联合先前得到的三弯矩方程组和两个补充方程得到恰定的三弯矩方程组。它是严格对角占优的三对角方程组，可以用追赶法求解。

结果：
a 直接计算即可得出
对于n=5有

对于n=10有

对于n=20有

b对于n=5有牛顿插值多项式

样条插值函数（分别对于各个小区间）为：

对于n=10有牛顿插值函数
p =

(9094987740384*X^9 + 5272349642029869237892763424*X^8 - 16938813146112*X^7 - 10950865478705002051580730624*X^6 + 10150693091136*X^5 + 7491941821973715378406714008*X^4 + 1915628554944*X^3 - 2014100801013926045821422321*X^2 + 322192441744*X + 210052147079480949741593257)/243810615467456022706126848

pp =

3.73*10^(-14)*X^9 + 21.6*X^8 - 6.95*10^(-14)*X^7 - 44.9*X^6 + 4.16*10^(-14)*X^5 + 30.7*X^4 + 7.86*10^(-15)*X^3 - 8.26*X^2 + 1.32*10^(-15)*X + 0.862

对于n=10有样条插值函数：

对于n=20有牛顿插值函数
- 6.82*10^(-10)*X^19 - 25770.0*X^18 + 2.52*10^(-9)*X^17 + 95660.0*X^16 - 3.79*10^(-9)*X^15 - 1.47*10^5*X^14 + 3.02*10^(-9)*X^13 + 1.21*10^5*X^12 - 1.46*10^(-9)*X^11 - 58660.0*X^10 + 4.3*10^(-10)*X^9 + 17220.0*X^8 - 6.6*10^(-11)*X^7 - 3066.0*X^6 + 4.06*10^(-12)*X^5 + 328.0*X^4 - 3.19*10^(-13)*X^3 - 21.6*X^2 + 3.08*10^(-15)*X + 0.993

对于n=20有样条插值函数：
- 5.25*X^3 - 14.37*X^2 - 12.98*X - 3.827
1.723*X^3 + 4.461*X^2 + 3.963*X + 1.256
0.07504*X^3 + 0.5066*X^2 + 0.7996*X + 0.4127
0.9738*X^3 + 2.394*X^2 + 2.121*X + 0.721
1.552*X^3 + 3.436*X^2 + 2.746*X + 0.846
3.552*X^3 + 6.434*X^2 + 4.245*X + 1.096
5.726*X^3 + 9.044*X^2 + 5.289*X + 1.235
12.54*X^3 + 15.17*X^2 + 7.127*X + 1.419
- 32.79*X^3 - 12.02*X^2 + 1.688*X + 1.056
- 89.07*X^3 - 28.91*X^2 - 2.345*10^(-15)*X + 1.0
89.07*X^3 - 28.91*X^2 - 6.395*10^(-15)*X + 1.0
32.79*X^3 - 12.02*X^2 - 1.688*X + 1.056
- 12.54*X^3 + 15.17*X^2 - 7.127*X + 1.419
- 5.726*X^3 + 9.044*X^2 - 5.289*X + 1.235
- 3.552*X^3 + 6.434*X^2 - 4.245*X + 1.096
- 1.552*X^3 + 3.436*X^2 - 2.746*X + 0.846
- 0.9738*X^3 + 2.394*X^2 - 2.121*X + 0.721
- 0.07504*X^3 + 0.5066*X^2 - 0.7996*X + 0.4127
- 1.723*X^3 + 4.461*X^2 - 3.963*X + 1.256
5.25*X^3 - 14.37*X^2 + 12.98*X - 3.827
c 容易计算各点函数值
YK =
0.0400 0.0416 0.0433 0.0451 0.0471 0.0491 0.0513 0.0536 0.0561 0.0588 0.0617 0.0648 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.8000 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.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 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.0152 -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 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 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.7749 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.7408 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.1839 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.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 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.1665
E (s)= 0.0031

这是由于已知数据的误差对插值多项式计算的结果带来了很大的误差。具体说，当函数值有误差而变为且误差限为。 由于插值与插值本质上是一样的，则插值与插值的误差估计是一样的。由插值多项式的计算值也会有误差，于是可以导出插值多项式的误差

当很大时，的误差就很大，于是就出现了前面观察到的现象。

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 + 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/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 - y6/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 - y4/22 - y5/36 - (5*y6)/78 - y7/28 - (7*y8)/90 - y9/48 - y2/54=y9

容易解得函数f（x）在各个节点处的值

用三次样条函数近似真解

90.779386156063992530107498168945*X^3 - 24.781023269490106031298637390137*X^2 - 0.00000000000363797880709171295166015625*X + 1.0000999999988380295690149068832

26.782958468349534086883068084717*X^3 + 19.304855964664966450072824954987*X^2 - 5.5107349042730220389785245060921*X + 1.2297139543440778197691543027759

5.7028477172279963269829750061035*X^3 - 5.0594986745181813603267073631287*X^2 + 0.58035375552140067156869918107986*X + 0.72212323269567946226743515580893

2.3772324005876726005226373672485*X^3 + 4.0305914580244461831171065568924*X^2 - 2.8284300441811751625209581106901*X + 1.1482212076578193205023126211017

0.18751811487527447752654552459717*X^3 + 0.74602002945584899862296879291534*X^2 - 1.1861443298974450044624973088503*X + 0.87450692194496326692387810908258

0.41370486008963780477643013000488*X^3 - 0.3812730486033615306951105594635*X^2 - 0.48158615610998367628781124949455*X + 0.72772396907193126480706268921494
-
3.2593013254809193313121795654297*X^3 + 7.8829908689303920255042612552643*X^2 - 6.6797840942605262171127833425999*X + 2.2772734536092258395001408644021

11.497100441840302664786577224731*X^3 - 30.852563770287815714254975318909*X^2 + 27.213826215055178181501105427742*X - 7.6083628866077219754515681415796

可以计算插值函数在节点处的值为
1.000100000000000 0.949494865095729 0.843069153460972
0.747070493078056 0.685497544329897 0.640100000000000
0.594851367010309 0.549968055375552 0.509600949631811
0.474997155522828 0.444500000000000 0.416600273048601
0.390883621796760 0.367218801767305 0.345690632989691
0.326600000000000 0.309798055670103 0.293299449779087
0.274890275110457 0.257659706332842 0.250000000000000

真解在节点处的值为

1.000000000000000 0.907029478458050 0.826446280991735
0.756143667296787 0.694444444444444 0.640000000000000
0.591715976331361 0.548696844993141 0.510204081632653
0.475624256837099 0.444444444444444 0.416233090530697
0.390625000000000 0.367309458218549 0.346020761245675
0.326530612244898 0.308641975308642 0.292184075967860
0.277008310249308 0.262984878369494 0.250000000000000

绿色为真解
红色为插值解可以看出他们非常相近

p229 第三题
用Newton迭代法求解方程组

算法分析：
用迭代法解非线性方程时，如何构造迭代函数是非常重要的，那么怎样构造的迭代函数才能保证迭代法收敛呢？牛顿迭代法就是常用的方法之一，其迭代格式的来源大概有以下几种方式：
1设，对在点作泰勒展开：

略去二次项，得到的线性近似式：。
由此得到方程0的近似根(假定0)，
即可构造出迭代格式(假定0)：
这就是牛顿迭代公式，若得到的序列{}收敛于，则就是非线性方程的根。
2 牛顿迭代法也称为牛顿切线法，这是由于的线性化近似函数＝是曲线＝过点的切线而得名的，求的零点代之以求的零点，即切线与轴交点的横坐标，如右图所示，这就是牛顿切线法的几何解释。实际上，牛顿迭代法也可以从几何意义上推出。利用牛顿迭代公式，由得到，从几何图形上看，就是过点作函数的切线，切线与轴的交点就是，所以有，整理后也能得出牛顿迭代公式：
。
3 要保证迭代法收敛，不管非线性方程0的形式如何，总可以构造：

作为方程求解的迭代函数。因为：
而且在根附近越小，其局部收敛速度越快，故可令：
若0(即根不是0的重根)，则由得：，
因此可令，则也可以得出迭代公式：。
4 迭代法的基本思想是将方程改写成等价的迭代形式，但随之而来的问题却是迭代公式不一定收敛，或者收敛的速度较慢。运用前述加速技巧，对于简单迭代过程，其加速公式具有形式：
，其中
记，上面两式可以合并写成：
这种迭代公式称作简单的牛顿公式，其相应的迭代函数是：。
需要注意的是，由于是的估计值，若取，则实际上便是的估计值。假设，则可以用代替上式中的，就可得到牛顿法的迭代公式：。
牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程逐步归结为某种线性方程来求解。对于方程组的解法和单个方程的解法类似。
第一个方程计算结果：

format long 格式后更精确地解如下

将所算结果代入方程可以得到误差满足程序设定的误差要求。
第二个方程的计算结果：

format long 格式后更精确地解如下

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

结果：

不妨减小步长，算出更为精确地结果 h=0.01时 结果如下图

function Gaussseidel
% 81页第五题用Gauss-seidel迭代求解电路电流方程组，使误差小于1e-3
clc
clear all
close all

A=[28 -3 0 0 0 ;-3 38 -10 0 -5; 0 -10 25 -15 0; 0 0 -15 45 0;0 -5 0 0 30]

b=[10 0 0 0 0]'
ep=1e-3;
x=maseidel(A,b,ep)
function x=maseidel(A,b,ep,x0)
% x0为初始迭代值，ep为最终解得误差
n=length(b);
if nargin<4;%若输入变量小于4个默认初始迭代值为零向量。
x0=zeros(n,1);
end
x=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;
end

function Jacobi
% 81页第五题用Jacobi迭代求解电路电流方程组，使误差小于1e-3
clear all
close all
clc

A=[28 -3 0 0 0 ;-3 38 -10 0 -5; 0 -10 25 -15 0; 0 0 -15 45 0;0 -5 0 0 30]

b=[10 0 0 0 0]'
ep=1e-3;
x=majacobi(A,b,ep)

function x=majacobi(A,b,ep,x0)
% x0为初始迭代值，ep为最终解得误差
n=length(b);
if nargin<4;% 若输入变量小于4个默认初始迭代值为零向量。
x0=zeros(n,1);
end
x=zeros(n,1);
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;
end

function zhuiganfa
clc
clear all
close 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)

function 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-