2022年西安交通大学-计算方法A2021上机实习报告 .pdf

1 XXXX学 院计算方法上机实验报告专业：班级：姓名：学号：日期：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 19 页2 一、用列主元 Gauss消去法解线性方程组Ax=b,其中：11111111821111111932111111114321111114,5432111118654321112376543211298765432136Ab1.1 算法组织消去法的中心就是“降维” ，即：将求解 n 元方程组的问题转化为先解n-1元方程组，一旦这个 n-1 元方程组的解取得，则剩余的一个未知量自然可以求得。这样逐步减少未知量个数的方法，便是求解多元方程组的一个重要思想。列主元消去法的基本思想是： 在进行第 k 步消元时，从第 k 列的对角线及其以下的各元素中选取绝对值最大的元素， 然后通过行变换将它交换到主元素的位置上，再进行消元。算法步骤如下：（1）选主元：在子块的第一列中选择一个元)( kkika使)(maxkiknikkkiaak并将第 k 行元与第ki行元互换。（2）消元计算：对 k=1,2,n-1 依次计算nkkibmbbnkkjiamaankkiaamkkikkikikkjikkijkijkkkkikkik,2, 1,2, 1,2, 1)()()1()()()1()()()(（3）回代求解)(1)()()()(iiinijjiijiiinnnnnnaxabxabx精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 19 页3 1 ,2,1 nni1.2 计算结果及分析运行 matlab程序，输入系数矩阵和行向量，运行结果如下图所示。求解结果和理论值相同，说明算法组织和matlab程序的正确性。二、求 20 阶三对角方程组 Tx=f 的解，其中：122021121(,),( 1,0,0, 1)12112TTxx xxTf2.1 算法组织系数矩阵T 是一种特殊的稀疏矩阵，在三次样条插值或者差分算法求解常微分方程边值问题中常会遇到，系数矩阵T 可以分解为一个上二对角阵和一个下二对角阵， 方程的求解问题转换为两个方程组的求解，相当于 Gauss消元法的部分称为“追”，相当于回代的过程称为“赶” ，追赶法的具体步骤如下：(1) 输入三对角矩阵T和右端向量f；(2) 预处理：将T压缩为四个一维数组iiiiabcf、，将分解矩阵压缩为三个一维数组iiiluy、(3) 追的过程：1111,ub yf精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 19 页4 2in111iiiiiiiiiiila uubl cyfl y(4) 回代求解，赶的过程nnxy1iiiikycxxu(5) 停止，输出结果2.2 计算结果及分析利用 matlab 编写追赶法程序，运行程序根据提示输入下对角元素、对角线元素、上对角元素以及右端的行向量，计算过程及结果如下图所示。经验证，计算结果和理论值相同，说明了算法和编程实现的准确性。三、用 Jacobi和 Gauss-Seidel方法求解线性方程组，某电流的电路方程满足方程组：12123523425342831033810501025150530015450iiiiiiiiiiiii试用 Jacobi和 Gauss-Seidel迭代法解之，使误差小于0.001。3.1 算法组织精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 19 页5 将原始线性代数方程组Axb改写为( )xx的形式，其中( )x为 x 的矩阵函数。于是可以得到迭代格式：(1)()()kkxx， 此即为 Jacobi迭代法的迭代格式。如果在计算(1)kx时，将已经算出的分量立即代换()kx对应分量，则得到Gauss Seidel迭代法的迭代格式。1. Jacobi迭代法的算法组织如下：(1) 给出迭代格式(1)()()kkxx(2) 给出迭代初始向量0 x、允许误差和最大迭代次数 N(3) 按照迭代格式(1)()()kkxx进行迭代，直至达满足迭代停止条件(1)()1maxkkiii nxx(4) 停止，输出结果2. Gauss Seidel迭代法的算法组织如下：(1) 给出迭代格式(1)()()kkxx(2) 给出迭代初始向量0 x、允许误差和最大迭代次数 N(3) 按照迭代格式(1)()()kkxx， 并且将已经算出的分量立即代换()kx对应分量进行迭代，直至达满足迭代停止条件(1)()1maxkkiii nxx(4) 停止，输出结果3.2 计算结果及分析利用 matlab 编写 Jacobi 迭代法程序，运行程序根据提示输入系数矩阵、右端行向量以及容许误差，计算过程及结果如下图所示。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 19 页6 利用 matlab编写 Gauss-Seidel迭代法程序，运行程序根据提示输入系数矩阵、右端行向量以及容许误差，计算过程及结果如下图所示。对比分析以上结果可知 ,达到相同的计算精度， Gauss Seidel 迭代比 Jacobi迭代的速度慢，迭代次数更多。四、Newton插值多项式和三次样条插值多项式已知21( )( 11)125f xxx，对5 , 10 , 20n1.计算函数( )f x在点21, (0,1, 2 ,)ixiinn处的值()if x; 2.求插值数据点,(0,1, 2, )iixyin的 Newton插值多项式( )nNx和三次样条插值多项式( )nSx; 3.对5 , 20n，计算21, (110,9099)100kxkk和相应的函数值,(),()knknkkyfxNxSx; 4.计算maxnnkkkE NyNx，maxnnkkkE SySx，解释所得到结果。4.1 算法组织4.1.1求 Newton插值多项式( )nNx，算法组织如下：Newton插值多项式的表达式如下：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 19 页7 010011( )()()()()nnnNxccxxcxxxxxx其中每一项的系数ci的表达式如下：12011010(,)(,)(,)iiiiifx xxf xxxcf xxxxx根据上述公式，为了得到系数需计算：(a) 一阶差商01,nf xf xf x(b) 二阶差商01121,nnf xxf xxf xx (c) n 阶差商01111,nnnf xxxf xxx(d) n+1阶差商011,nnf xxxx4.1.2求三次样条插值多项式( )nSx，算法组织如下：所谓三次样条插值多项式( )nSx是一种对区间进行分段的分段函数，然后在每一段上进行分析，即它在节点ix011()nnaxxxxb分成的每个小区间1,iixx上是 3 次多项式，其在此区间上的表达式如下：22331111111( )()()()()666,1,2, .iiiiiiiiiiiiiiiiihxxhxxS xxxMxxMyMyMhhhxxxin，因此，只要确定了iM的值，就确定了整个表达式，iM的计算方法如下：令：11111111116()6 (,)iiiiiiiiiiiiiiiiiiiiihhhhhhyyyydfxx xhhhh，则iM满足如下 n-1 个方程：1121,2,1iiiiiiMMMd in，方程中有 n+1个未知量， 则令0M和nM分别为零，则由上面的方程组可得到精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 19 页8 (11)iMin的值，可得到整个区间上的三次样条插值多项式( )nSx。4.2 计算结果及分析第一部分运行结果第二部分运行结果：n=5时，Newton 插值多项式Nnn=- 0.000000000000003664*xx5 + 1.202*xx4 + 0.000000000000001465*xx3 - 1.731*xx2 - 0.00000000000000005276*xx + 0.5673 n=10时，Newton 插值多项式Nnn=- 220.9*xx10 - 0.00000000000005684*xx9 + 494.9*xx8 + 0.0000000000002187*xx7 - 381.4*xx6 + 0.0000000000000964*xx5 + 123.4*xx4 - 0.00000000000002632*xx3 - 16.86*xx2 + 0.0000000000000008847*xx + 1.0 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 19 页9 n=20时，Newton 插值多项式Nnn=260220.0*xx20 - 1012200.0*xx18 + 0.00000000002561*xx17 + 1639900.0*xx16 - 0.000000000005128*xx15 - 1443300.0*xx14 + 0.000000000007088*xx13 + 757330.0*xx12 - 0.00000000008592*xx11 - 245220.0*xx10 + 0.00000000001925*xx9 + 49322.0*xx8 + 0.000000000003205*xx7 - 6119.0*xx6 + 0.0000000000001079*xx5 + 470.8*xx4 - 0.0000000000001146*xx3 - 24.14*xx2 + 0.00000000000000572*xx + 1.0 n=5时，三次样条插值多项式f = - 1.965*t3 - 5.06504*t2 - 4.16107*t - 1.02257, 8.57202*t3 + 16.009*t2 + 9.88828*t + 2.09951, - 21.1146*t3 - 13.6776*t2 - 0.00726356*t + 1.0, 21.2454*t3 - 13.6776*t2 - 0.00726356*t + 1.0, - 9.22574*t3 + 16.7935*t2 - 10.1643*t + 2.12856 n=10时，三次样条插值多项式f = - 0.0694247*t3 - 0.0604978*t2 + 0.161278*t + 0.190813, 1.3153*t3 + 3.33838*t2 + 2.94218*t + 0.94924, - 1.33175*t3 - 1.71509*t2 - 0.273666*t + 0.267091, 16.726*t3 + 22.9092*t2 + 10.9192*t + 1.96298, - 36.235*t3 - 20.4226*t2 - 0.898565*t + 0.888637, 0.00534817*t3 - 10.5388*t2 - 0.0000441997*t + 0.915865, 36.2136*t3 - 20.4138*t2 + 0.897681*t + 0.888661, - 16.6512*t3 + 22.8392*t2 - 10.8986*t + 1.96105, 1.05365*t3 - 1.30377*t2 + 0.0754743*t + 0.298313, - 0.277758*t3 + 1.23801*t2 - 1.54202*t + 0.641417 n=20时，三次样条插值多项式f = 0.166208*t3 + 0.601713*t2 + 0.778802*t + 0.381759, 0.260135*t3 + 0.856657*t2 + 1.00947*t + 0.451324, 0.439278*t3 + 1.29172*t2 + 1.36166*t + 0.54636, 0.766898*t3 + 1.99376*t2 + 1.86312*t + 0.665756, 1.43578*t3 + 3.23598*t2 + 2.63211*t + 0.824436, 2.70965*t3 + 5.23777*t2 + 3.68067*t + 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 19 页10 1.00752, 5.55772*t3 + 8.89958*t2 + 5.25001*t + 1.23171, 8.44392*t3 + 11.7858*t2 + 6.21208*t + 1.33861, 2.75729*t3 + 7.72389*t2 + 5.24496*t + 1.26185, - 93.2743*t3 - 33.4325*t2 - 0.634519*t + 0.981876, - 0.0000138825*t3 - 20.1076*t2 + 0.0000000314796*t + 0.991947, 93.2743*t3 - 33.4325*t2 + 0.634519*t + 0.981876, - 2.75748*t3 + 7.72399*t2 - 5.24498*t + 1.26185, - 8.4432*t3 + 11.7852*t2 - 6.21194*t + 1.33859, - 5.56042*t3 + 8.90243*t2 - 5.25101*t + 1.23182, - 2.6996*t3 + 5.22424*t2 - 3.67464*t + 1.00663, - 1.47329*t3 + 3.29719*t2 - 2.66523*t + 0.830383, - 0.626907*t3 + 1.72533*t2 - 1.69218*t + 0.629594, - 0.961731*t3 + 2.44281*t2 - 2.20466*t + 0.751614, 1.68969*t3 - 3.99635*t2 + 3.00799*t - 0.654974 第三部分运行结果：n=5时，从左至右，函数值，牛顿插值，三次样条插值n=20时，从左至右，函数值，牛顿插值，三次样条插值精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 19 页11 第四部分运行结果当 n=5 时：ENn = 0.432692307692308 ；ESn = 0.128807360154696 ；当 n=10 时：ENn =1.915643050219580 ；ESn = 0.084135000000000 ；当 n=20 时：ENn = 58.278125107950132 ；ESn = 0.045126353790614 ；由上面的结果显示，相同的n 值下， Newton 插值多项式误差比三次样条多项式要打的多；使用Newton 插值多项式出现随着n 的增大，误差也逐渐增大的现象，其最大误差达到58.27813，而相应的三次样条插值多项式随着n 的增大，误差逐渐减小， n=20 的误差仅为 0.003088 。可见， n 越大， Newton 插值越可能偏离被插值函数， 而相应的三次样条插值则能更接近于被插值函数。因此，三次样条更合适进行插值计算。五、分别用 n=2，4，8，16 的复化 Simpson公式求积分方程3212045251()( )( )8 (1)1xxxxy tdty xxt的近似解。然后用三次样条插值函数近似y(t)的离散近似解，并在区间0,1的一精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 19 页12 些点上与真解21( )(1)y xx比较。5.1 算法组织复化梯形公式的递推算法求积分的主要思想是利用若干小梯形的面积代替原方程的积分，当精度达不到要求时， 可以通过增加点数对已有的区间再次划分，达到所需精度时即可。5.2 计算结果及分析梯形求积 n = 2时求解的离散值为：0.000000 1.098016 0.500000 0.498214 1.000000 0.259524 梯形求积 n = 4时求解的离散值为：0.000000 1.024813 0.250000 0.659179 0.500000 0.457989 0.750000 0.334440 1.000000 0.252275 梯形求积 n = 8时求解的离散值为：0.000000 1.006236 0.125000 0.795651 0.250000 0.644819 0.375000 0.533035 0.500000 0.447845 0.625000 0.381390 0.750000 0.328513 0.875000 0.285718 1.000000 0.250565 0.000000 1.001561 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 19 页13 0.062500 0.887286 0.125000 0.791507 0.187500 0.710436 0.250000 0.641206 0.312500 0.581616 0.375000 0.529954 0.437500 0.484872 0.500000 0.445296 0.562500 0.410362 0.625000 0.379372 0.687500 0.351751 0.750000 0.327027 0.812500 0.304807 0.875000 0.284763 0.937500 0.266619 1.000000 0.250141 和三次样条结果对比分析，如下图所示：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 19 页14 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 19 页15 从上面的结果图可以看出， 随着 n 的增大，梯形求积分结果准确度增加，当n=16时，结果已经非常准确。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 19 页16 六、用 Romberg 方法计算定积分：0.830.3sin( )xxIdxx，(0)1f使误差不超过510。6.1 算法组织在变步长积分的基础上，可以证明变步长积分的结果221()3nnnIfTTT与复化Simpson求积公式的关系为：221()41nnnnSIfTTT同理可以证明222141nnnnCSSS223141nnnnRCCC如此进行下去即为Romberg积分方法，以上称为Romberg递推公式。6.2 计算结果及分析利用 matlab 编写程序，运行程序根据提示输入参数，计算过程及结果如下图所示。七、用 Newton迭代法解如下的方程组22222T2221.0241.0 ,341.0 xyzxyzxyz初始向量为（1.0 ，1.0 ，1.0 ）7.1 算法组织(a) 先构造牛顿迭代函数，( ) /( )rxf xfx。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 16 页，共 19 页17 (b) 设置初始值(0)x和精度范围 eps。(c) 进行第一次迭代(0)(0)(0)() /()rxf xfx。(d) 判断(0)|rxeps是否成立，若成立则迭代完成， 若不成立则继续迭代。(e) 迭代直到(0)|rxeps成立迭代完成。7.2 计算结果及分析利用 matlab 编写程序，运行程序根据提示输入参数，计算过程及结果如下图所示。结果分析，计算结果在误差范围内，符合要求。八、用 Runge-Kutta 方法计算高阶方程的初值问题。23(0)1,(0)3,(0)2yyyytyyy分别去步长 h=1/2，1/8, 1/32，1/128,计算 y(1)并与解析解( )21ty ttet做比较。8.1 算法组织1.输入待求微分方程( ,( )f t y t，自变量的范围( , )a b，初值(0)y，步长 h，以及求解范围内的取点数 n。2.求解精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 17 页，共 19 页18 1123412132431(22)6(,)1(,)221(,)22(,)iiiiiiiiiiyyKKKKKhf x yhKhf xyKhKhf xyKKhf xh yK8.1 计算结果及分析利用 matlab 编写程序，运行程序根据提示输入参数，计算过程及结果如下图所示。从上边计算结果可以看出，随着步长的减小，计算值和理论值误差快速降低。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 18 页，共 19 页19 九、总结通过对本次上机题，学习和掌握了基本的matlab 编程方法和技巧，以及matlab在数值求解和计算中的优点； 在完成各个题目的过程中， 深化了理论知识的学习，编写 matlab 算法实现了运算过程的自动化，并不通过不同参数变化对结果影响的对比分析，进一步认识了不同算法的实际应用。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 19 页，共 19 页