中 中 南 南 林 林 业 业 科 科 技 技 大 大 学 学 电 电 子 子 与 与 信 信 息 息 工 工 程 程 学 学 院 院 《数值分析》 实 实 验 验 指 指 导 导 书 书 陈爱斌 编 计算机教研室 2006年2月 目 录 前 言................................................... 1 实验一 函数插值方法...................................... 2 实验二 函数逼近与曲线拟合................................ 3 实验三 数值积分与数值微分................................ 4 实验四 线方程组的直接解法................................ 5 实验五 解线性方程组的迭代法.............................. 7 实验六 非线性方程求根.................................... 8 实验七 矩阵特征值问题计算................................ 9 实验八 常微分方程初值问题数值解法....................... 11 附录一：实验报告内容要求................................ 12 附录二：部分程序示例.................................... 13 1、〖多项式系列程序〗................................. 13 2、〖最小二乘法曲线拟合程序〗......................... 15 3、〖Romberg算法计算数值积分程序〗................... 18 4、〖线方程组的直接解法系列程序〗..................... 19 5、〖Jacobi迭代法解线性方程组程序〗.................. 24 6、〖牛顿法非线性方程求根程序〗....................... 25 7、〖幂法求矩阵特征值程序〗........................... 26 8、〖四阶经典的龙格-库塔方法解常微分方程程序〗........ 27 数值分析实验指导书 1 前 言 第1页 结合课程教学，配备适当的上机实验(16个学时)以便加深课堂 教学的实践性，同时通过实验可以加强对数学模型的总体分析，算法 选取，程序结构，上机调试和结果分析等环节的训练。本实验指导书 共包含8个实验，要求学生在16个实验课时内完成。为使实验更为 有成效,需要写出实验报告（格式要求见附录）,以此可作为《数值分 析》课程成绩评定的参考。 数值分析实验指导书 2 实验一 函数插值方法 一、问题提出 对于给定的一元函数 ( )xfy = 的n+1个节点值 ( ) jj xf y = 。 试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。 ()nj ,,1,0L= 数据如下： （1） j x 0.4 0.55 0.65 0.80 0.95 1.05 j y 0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次Lagrange多项式L，和分段三次插值多项式，计算()x 5 ( )(99.0,596.0 ff ) 的 值。（提示：结果为 ()≈596.0f 0.625732 ( )≈99.0f 1.05423 ） （2） j x 1 2 3 4 5 6 7 j y 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式L 6 ，计算( )8.1f的值。 ( )x 结果()≈8.1f 0.164762 ( )≈15.6f 0.001266 二、要求 1、 利用Lagrange插值公式 () k n ki i ik i n k n y xx xx xL ? ? ? ? ? ? ? ? ? ? ? ? ∑= ∏ ≠ = = 0 0 编写出插值多项式程序； 2、 给出插值多项式或分段三次插值多项式的表达式； 3、 根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何； 4、 对此插值问题用Newton插值多项式其结果如何。Newton插值多项式如下： () ∏ ? ≠ = = ??∑+= 1 0 0 1 0 )(],[)( k kj j jk n k n xxxxfxfxNL 其中： ∏ ≠ = = ? ∑= k ij j ji i k i k xx xf xxf 0 0 0 )( )( ],,[L 三、目的和意义 1、 学会常用的插值方法，求函数的近似表达式，以解决其它实际问题； 2、 明确插值多项式和分段插值多项式各自的优缺点； 3、 熟悉插值方法的程序编制； 第2页 4、 如果绘出插值函数的曲线，观察其光滑性。 数值分析实验指导书 3 实验二 函数逼近与曲线拟合 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验 中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟 合曲线。 y 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为； () 3 3 2 21 tatatat ++=? 3、打印出拟合函数()t?，并打印出( ) j t?与( ) j ty的误差，12,,2,1L=j； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、* 绘制出曲线拟合图。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系 第3页 (t分） 0 5 10 15 20 25 30 35 40 45 50 55 ( ) 4 10 ? ×y 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 数值分析实验指导书 4 实验三 数值积分与数值微分 一、问题提出 选用复合梯形公式，复合Simpson公式，Romberg算法，计算 （1） I = dxx ∫ ? 4 1 0 2 sin4 ( )5343916.1≈I （2） I = dx x x ∫ 1 0 sin ()9460831.0,1)0( ≈= I f （3） I = dx x e x ∫ + 1 0 2 4 （4） I = () dx x x ∫ + + 1 0 2 1 1ln 二、要求 1、 编制数值积分算法的程序； 2、 分别用两种算法计算同一个积分，并比较其结果； 3、 分别取不同步长( / ab h )?= n，试比较计算结果（如n = 10, 20等）； 4、 给定精度要求ε，试用变步长算法，确定最佳步长。 三、目的和意义 1、 深刻认识数值积分法的意义； 2、 明确数值积分精度与步长的关系； 第4页 3、 根据定积分的计算方法，可以考虑二重积分的计算问题。 数值分析实验指导书 5 实验四 线方程组的直接解法 一、问题提出 给出下列几个不同类型的线性方程组，请用适当算法计算其解。 1、设线性方程组 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ????? ? ??? ?? ?? ??? ??? ??? ?? ?? 13682438100 412029137264 221234179111016 1035243120 53621775868 3233761624 4911315120 1301231224 0010563568 0000121324 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 10 9 8 7 6 5 4 3 2 1 x x x x x x x x x x = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 21 19 38 13 46 3 2 3 12 5 x = ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 ) T * 2、设对称正定阵系数阵线方程组 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ?? ??? ???? ???? ???? ?? ? 192433600 2141103520 411144334 3104221812 33416120 653811414 02312122 00420424 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 8 7 6 5 4 3 2 1 x x x x x x x x = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 45 15 22 9 23 20 6 0 x = ( 1, -1, 0, 2, 1, -1, 0, 2 ) T * 第5页 数值分析实验指导书 6 3、三对角形线性方程组 = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ?? ?? ?? ?? ?? ?? ?? ? 4100000000 1410000000 0141000000 0014100000 0001410000 0000141000 0000014100 0000001410 0000000141 0000000014 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 10 9 8 7 6 5 4 3 2 1 x x x x x x x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 5 5 4 14 12 6 2 13 5 7 x = ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 ) T * 二、要求 1、 对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与 改进平方根法；追赶法求解（选择其一）； 2、 应用结构程序设计编出通用程序； 3、 比较计算结果，分析数值解误差的原因； 4、 尽可能利用相应模块输出系数矩阵的三角分解式。 三、目的和意义 1、通过该课题的实验，体会模块化结构程序设计方法的优点； 2、运用所学的计算方法，解决各类线性方程组的直接算法； 3、提高分析和解决问题的能力，做到学以致用； 第6页 4、通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。 数值分析实验指导书 7 实验五 解线性方程组的迭代法 一、问题提出 对实验四所列目的和意义的线性方程组，试分别选用Jacobi 迭代法，Gauss-Seidol 迭代法和SOR方法计算其解。 二、要求 1、体会迭代法求解线性方程组，并能与消去法做以比较； 2、分别对不同精度要求，如由迭代次数体会该迭代法的收敛快慢； 543 10,10,10 ??? =ε 3、对方程组2，3使用SOR方法时，选取松弛因子ω =0.8，0.9，1，1.1，1.2等，试 看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者； 4、给出各种算法的设计程序和计算结果。 三、目的和意义 1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较； 2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序； 3、体会上机计算时，终止步骤 ∞ + ? )()1( kk xx < ε 或k >（予给的迭代次数），对迭 代法敛散性的意义； 第7页 4、体会初始解 x ( ，松弛因子的选取，对计算结果的影响。 )0 数值分析实验指导书 8 实验六 非线性方程求根 设方程f(x)=x - 3x –1=0 有三个实根 x * 1 =1.8793 , x * =-0.34727 ,x =-1.53209 现采用下面六种不同计算格式，求 f(x)=0的根 x * 1 或x * 2 3 2 * 3 一、问题提出 1、 x = 2 13 x x + 2、 x = 3 1 3 ?x 3、 x = 3 13 +x 4、 x = 3 1 2 ?x 5、 x = x 1 3+ 6、 x = x - ( ) 1 13 3 1 2 3 ? ?? x xx 二、要求 1、编制一个程序进行运算，最后打印出每种迭代格式的敛散情况； 2、用事后误差估计 kk xx ? +1 ? ε来控制迭代次数，并且打印出迭代的次数； 3、初始值的选取对迭代收敛有何影响； 4、分析迭代收敛和发散的原因。 三、目的和意义 1、通过实验进一步了解方程求根的算法； 2、认识选择计算格式的重要性； 3、掌握迭代算法和精度控制； 第8页 4、明确迭代收敛性与初值选取的关系。 数值分析实验指导书 9 实验七 矩阵特征值问题计算 一、问题提出 利用冪法或反冪法，求方阵A=(a ) n 的按模最大或按模最小特征值及其对应的特 征向量。 ij n× 设矩阵A的特征分布为： nn λλλλλ ≥≥≥? ?1321 L 且Ax = j jj xλ 试求下列矩阵之一 （1） A= 求 ? ? ? ? ? ? ? ? ? ? ? ? ? 611 142 121 1 λ，及x 1 () 取= ( 1, 1, 1 ) , 0 υ T ε = 10 5? 结果 λ 1 ≈ -6.42106, x 1 ≈ ( -0.046152, -0.374908, 1 ) T （2） A= 求 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? 421578 235341 156213 532717 741152 813724 1 λ , 6 λ 及 x 1 () 取 0 υ T )1,0,0,1,0,1(≈ 5 10 ? =ε 结果： ( ) T 16 1.0 0.4972, 0.5644, 0.9973, 0.5401, 0.8724, x 1.62139 ≈≈≈ λλ 30525.21 1 （3） A= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ?? ?? ? 21 121 121 121 12 求 1 λ 及 x1 取 ( )0 υ =( 1, 1, 1, 1, 1 ) T 4 10 ? =ε 结果≈λ 3.7321 第9页 数值分析实验指导书 10 （4） A= ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? 1254 2613 5131 4312 取 ( )0 υ =( 1, 1, 1, 1 ) T 2 10 ? =ε 这是一个收敛很慢的例子，迭代1200次才达到10 5? 结果 1 λ ≈ -8.02857835 x 1 T )564212.2,757730.0,501460.2,1( ??≈ （5） A= ? ? ? ? ? ? ? ? ? ? ? ? ? 611 142 121 有一个近似特征值 –6.42，试用反冪法求对应的特征向量，并改进特征值（原点 平移法）。 取 ()0 υ = ( 1, 1, 1 ) T 4 10 ? =ε 结果 ≈λ -6.42107 x T )1,37918.0,0461465.0( ??≈ 二、要求 1、掌握冪法或反冪法求矩阵部分特征值的算法与程序设计； 2、会用原点平移法改进算法，加速收敛；对矩阵B=A-PI取不同的P值，试求其效 果； 3、试取不同的初始向量，观察对结果的影响； ()0 υ 4、对矩阵特征值的其它分布，如 21 λλ = 且 321 λλλ ≥=如何计算。 三、目的和意义 1、求矩阵的部分特征值问题具有重要实际意义，如求矩阵谱半径()Aρ =max 1 λ， 稳定性问题往往归于求矩阵按模最小特征值； 2、进一步掌握冪法、反冪法及原点平移加速法的程序设计技巧； 第10页 3、问题中的题（5），反应了利用原点平移的反冪法可求矩阵的任何特征值及其特征 向量。 数值分析实验指导书 11 实验八 常微分方程初值问题数值解法 一、问题提出 科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler 法，Rung-Kutta方法求其数值解，诸如以下问题： （1） () ? ? ? ? = ?=′ 0 y xy y x y 0 4 ? 0 < x≤2 分别取h=0.1,0.2,0.4时数值解。 初值问题的精确解 2 54 x ey ? +=。 （2）用r=3的Adams显式和预 - 校式求解 () ? ? ? ? ? ?=′ =? 22 01 yxy y 01 ≤≤? x 取步长h=0.1，用四阶标准R-K方法求值。 （3）用改进Euler法或四阶标准R-K方法求解 ( ) () () ? ? ? ? ? =?=′ =?=′ ?=′=′ 10 33 00 12 10 121 3 y yy 2 y yy y y y 10 ≤≤ x 取步长0.01,计算=h ( ) ( ) ( )15.0,10.0,05.0 y y y数值解，参考结果 () ( ) ( ) 8613125.015.0,9880787.015.0 1 ≈,1493359.015.0 32 ≈?≈ y y y （4）利用四阶标准R- K方法求二阶方程初值问题的数值解 （I） () () ? ? ? =′= =+′?′′ 10,00 023 yy yyy 0.02h x =≤≤ 10 (II) () () ? ? ? =′= =+′??′′ 00,10 0)1(1.0 2 y y yyyy 0.1h x =≤≤ ,10 (III) () () ? ? ? ? ? =′= + =′ 00,10 1 yy e y y x 0.1h x =≤≤ ,20 (IV) ? 0 () () ? ? =′= =+′′ 00,10 0sin yy yy 2,4 0.h x =≤≤ 二、要求 1、 根据初值问题数值算法，分别选择二个初值问题编程计算； 2、 试分别取不同步长，考察某节点处数值解的误差变化情况； j x 3、 试用不同算法求解某初值问题，结果有何异常； 4、 分析各个算法的优缺点。 三、目的和意义 1、 熟悉各种初值问题的算法，编出算法程序； 2、 明确各种算法的精度与所选步长有密切关系； 第11页 3、 通过计算更加了解各种算法的优越性。 数值分析实验指导书 12 附录一：实验报告内容要求 实验报告内容要求 一、课题名称 二、班级、姓名 三、目的和意义 方法的理论意义和实用价值，如解超越方程的弦截法，改进 了牛顿法，它适用于任意连续函数在大范围中求解，并且避免计 算导数值，使其更具有实用性。 四、计算公式 五、结构程序设计 六、结果讨论和分析 第12页 如初值对结果的影响；不同方法的比较；该方法的特点和改 进；整个实验过程中（包括程序编写，上机调试等）出现的问题 及其处理等广泛的问题，以此扩大知识面和对实验环节的认识。 数值分析实验指导书 13 附录二：部分程序示例 （本部分程序与本指导书配套而且均已采用Turbo C 2.0调试通过） 1、〖多项式系列程序〗 (1) Lagrange插值 第13页 #define N 5 float x[] = {0.4, 0.55, 0.65, 0.80, 0.95, 1.05}; float y[] = {0.41075, 0.57815, 0.69675, 0.90, 1.00, 1.25382}; float p(float xx) { int i,k; float pp = 0, m1, m2; for( i=0; i<=N; i++ ) { m1 = 1; m2 = 1; for( k=0; k<=N; k++ ) if( k != i ) { m1 *= xx - x[k]; m2 *= x[i] - x[k]; } pp += y[i] * m1/m2; } return pp; } main() { printf( "f(0.596)=%lf\n", p(0.596) ); printf( "f(0.99)=%lf\n", p(0.99) ); } 数值分析实验指导书 14 (2) 牛顿插值 第14页 #define N 6 float x[]={1,2,3,4,5,6,7}; float y[]={0.368,0.135,0.050,0.018,0.007,0.002,0.001}; float jc( int k ) /* 此函数实现求均差 f[x 0 ,x 1 ,...,x k ] */ { int i, j; float s = 0, m; for(i=0; i<=k; i++) { for(m=1,j=0; j<=k; j++) if( j!=i ) m *= x[i]-x[j]; s += y[i] / m; } return s; } float newton( float X ) /* 此函数实现求牛顿插值函数N n (x) */ { float nt=y[0], m; int i, j; for( i=1; i<=N; i++) { for( m=1,j=0; j<=i-1; j++) m *= X-x[j]; nt += jc(i) * m; } return nt; } main() { float X=1.8; printf("%lf",newton(X)); } 数值分析实验指导书 15 2、〖最小二乘法曲线拟合程序〗 第15页 /*曲线拟合*/ #include "stdio.h" #include "math.h" #define num 10 float neiji(float b[num],float c[num]) /*内积函数*/ { int p; float nj=0; for (p=1;p<num;p++) nj+=c[p]*b[p]; return nj; } float s[num],x[num],y[num],fai[num][num],afa[num]; float beida[num],a[num],xfai[num],yd[num],max,pcpfh; void main() { int i,j,k,n,index,flag; char conti; conti=' '; printf("请输入已知点的个数n=\n"); scanf("%d",&n); printf("请输入x和y:"); for(i=1;i<=n;i++) { printf("x[%d]=",i); scanf("%f",&x[i]); printf("y[%d]=",i); scanf("%f",&y[i]); } 数值分析实验指导书 16 第16页 while(conti==' ') { printf("请输入拟和次数="); scanf("%d",&index); pcpfh=0; afa[1]=0; a[0]=0; for(i=1;i<=n;i++) {afa[1]+=x[i]; a[0]+=y[i]; fai[0][i]=1; } afa[1]=afa[1]/n; a[0]=a[0]/n; for (i=1;i<=n;i++) { fai[1][i]=x[i]-afa[1]; } a[1]=neiji(fai[1],y)/neiji(fai[1],fai[1]); for(k=1;k<index;k++) { for(i=1;i<=n;i++) xfai[i]=x[i]*fai[k][i]; afa[k+1]=neiji(fai[k],xfai)/neiji(fai[k],fai[k]); beida[k]=neiji(fai[k],fai[k])/neiji(fai[k-1],fai[k-1]); for(j=1;j<=n;j++) fai[k+1][j]=(x[j]-afa[k+1])*fai[k][j]-beida[k]*fai[k-1][j]; a[k+1]=neiji(fai[k+1],y)/neiji(fai[k+1],fai[k+1]); } printf("%d次拟和结果为\n",index); for(i=0;i<=index;i++) printf("a[%d]=%f\n",i,a[i]); for(i=1;i<=index;i++) printf("afa[%d]=%f\n",i,afa[i]); for(i=1;i<index;i++) printf("beida[%d]=%f\n",i,beida[i]); 数值分析实验指导书 17 第17页 < 接上页 > for(i=1;i<=n;i++) { for(k=0;k<=index;k++) s[i]+=a[k]*fai[k][i]; yd[i]=(float)(fabs(y[i]-s[i])); pcpfh+=yd[i]*yd[i]; s[i]=0; } max=0; for(i=1;i<=n;i++) if(yd[i]>max) {max=yd[i]; flag=i; } printf("当x=%f时,偏差最大=%f,偏差平方和为%f\n",x[flag],max,pcpfh); printf("继续拟和请按space,按其他键退出"); conti=getchar(); conti=getchar(); } } 数值分析实验指导书 18 3、〖Romberg算法计算数值积分程序〗 第18页 float f(float x) { return( exp( x ) / ( 4 + x*x ) ); } main() { float a=0,b=1; float h=b-a,T1,T2,s,x; int i; T1=h/2*(f(a)+f(b)); for(i=0;i<3;i++) { printf("h=%f,T=%f\n",h,T1); s=0; x=a+h/2; while(x<b) { s+=f(x); x+=h; } T2=T1/2+h/2*s; h/=2; T1=T2; } printf("h=%f,T=%f\n",h,T1); return; } 数值分析实验指导书 19 4、〖线方程组的直接解法系列程序〗 以下分别用Gauss消去法、平方根法、追赶法编程解线性方程组。这 几个程序中，数据文件格式分别如下： 8 4 2 2 -4 -1 14 0 -2 1 6 2 1 -8 - 4 3 -3 - 0 2 5 - 0 0 6 3 sqrt.dat 第19页 10 4 -1 -1 4 -1 4 -1 4 -1 4 -1 4 -1 4 -1 4 -1 4 -1 4 4 2 -3 -1 2 1 0 0 0 0 5 8 6 -5 -3 6 5 0 1 0 0 1 4 2 -2 -1 3 2 -1 0 3 1 3 0 -2 1 5 -1 3 -1 1 9 4 2 -4 2 6 -1 6 7 -3 3 2 3 3 8 6 -1 3 7 17 2 6 -3 5 4 0 2 -1 3 -4 2 5 3 0 1 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 - 1 22 4 4 11 3 -10 1 14 -3 -4 2 19 zgf.dat Gauss.dat 2 6 3 38 19 21 7 -1 5 -1 -13 -1 2 -1 6 -1 -12 -1 14 -1 -4 -1 5 -5 数值分析实验指导书 20 (1) Gauss消去法解线性方程组 第20页 #include "stdio.h" main() { FILE *f; double a[15][15]; double b[15],s; int i,j,k,n; f=fopen("Gauss.dat","r"); fscanf(f,"%d",&n); for(i=1;i<=n;i++) { for(j=1;j<=n;j++) fscanf(f,"%lf",&a[i][j]); fscanf(f,"%lf",&b[i]); } fclose(f); k=1; do { for (j=k+1;j<=n;j++) a[k][j]=a[k][j]/a[k][k]; b[k]=b[k]/a[k][k]; i=1; for(i=k+1;i<=n;i++) { for (j=k+1;j<=n;j++) a[i][j]=a[i][j]-a[i][k]*a[k][j]; b[i]=b[i]-a[i][k]*b[k]; } if(k==n) break; k++; }while(1); for(i=n-1;i>=1;i--) { s=0; for(j=i+1;j<=n;j++) s=s+a[i][j]*b[j]; b[i]=b[i]-s; } for(i=1;i<=n;i++) printf("b[%2d]=%lf\n",i,b[i]); } 数值分析实验指导书 21 (2)平方根法求三角阵 第21页 /*************A=L*L T ***************/ #include "stdio.h" #include "math.h" main() { FILE *f; double a[15][15],l[15][15]; double s; int i,j,k,n; f=fopen("sqrt.dat","r"); fscanf(f,"%d",&n); for(i=1;i<=n;i++) for(j=1;j<=i;j++) { fscanf(f,"%lf",&a[i][j]); } fclose(f); for(i=1;i<=n;i++) { for(j=1;j<=i-1;j++) { for(s=0,k=1;k<=j-1;k++) s+=l[i][k]*l[k][j]; l[i][j]=(a[i][j]-s)/l[j][j]; } s=0; for(s=0,k=1; k<=i-1; k++) s+=l[i][k]*l[i][k]; l[i][i]=sqrt(a[i][i]-s); } printf("\n=======================\n"); for(i=1;i<=n;i++) { for(j=1;j<=i;j++) printf("%lf ", l[i][j]); printf("\n"); } } 数值分析实验指导书 22 (3)改进的平方根法求解线性方程组 第22页 #include "stdio.h" main() { FILE *f; double a[15][15],b[15],d[15],l[15][15],x[15],y[15]; double s; int i,j,k,n; /******************************************************/ f=fopen("sqrt.dat","r"); fscanf(f,"%d",&n); for(i=1;i<=n;i++) { for(j=1;j<=i;j++) fscanf(f,"%lf",&a[i][j]); fscanf(f,"%lf",&b[i]); } fclose(f); /*****************************************************/ for(i=1;i<=n;i++) { for(j=1;j<=i-1;j++) { for(s=0,k=1;k<=j-1;k++) s+=d[k]*l[i][k]*l[j][k]; l[i][j]=(a[i][j]-s)/d[j]; } s=0; for(s=0,k=1; k<=i-1; k++) s+=d[k]*l[i][k]*l[i][k]; d[i]=a[i][i]-s; } /***************************************************/ for(i=1;i<=n;i++) { for(s=0,k=1;k<=i-1;k++) s+=l[i][k]*y[k]; y[i]=b[i]-s; } /****************************************************/ for(i=n;i>=1;i--) { for(s=0,k=i+1;k<=n;k++) s+=l[k][i]*x[k]; x[i]=y[i]/d[i]-s; } /*****************************************************/ for(i=1;i<=n;i++) printf("x[%d]=%lf\n",i,x[i]); } 数值分析实验指导书 23 (4)追赶法求解线性方程组 第23页 #include "stdio.h" main() { FILE *f; double a[15],b[15],c[15],d[15]; double t; int i,n; /**********************************************/ f=fopen("zgf.dat","r"); fscanf(f,"%d",&n); fscanf(f,"%lf%lf%lf",&b[1],&c[1],&d[1]); for(i=2;i<=n-1;i++) { fscanf(f,"%lf%lf%lf%lf",&a[i],&b[i],&c[i],&d[i]); } fscanf(f,"%lf%lf%lf",&a[n],&b[n],&d[n]); fclose(f); /*********************************************/ c[1]=c[1]/b[1]; d[1]=d[1]/b[1]; for(i=2;i<=n-1;i++) { t=b[i]-c[i-1]*a[i]; c[i]=c[i]/t; d[i]=(d[i]-d[i-1]*a[i])/t; } d[n]=(d[n]-d[n-1]*a[n])/(b[n]-c[n-1]*a[n]); for(i=n-1;i>=1;i--) d[i]=d[i]-c[i]*d[i+1]; printf("\n********************************\n"); for(i=1;i<=n;i++) printf("d[%2d]=%lf\n",i,d[i]); } 数值分析实验指导书 24 5、〖Jacobi迭代法解线性方程组程序〗 第24页 #define M 3 #include "math.h" main() { double a[M][M]={ {10,-1,-2}, {-1,10,-2}, {-1,-1, 5} }; double b[M] = { 7.2, 8.3, 4.2 }; double x[M] = {0,0,0},y[M]; double s, max, eps = 0.000000001; int k, i, j, N = 100; k = 1; while( k<N ) { for ( i=0; i<M ; i++ ) { s = 0; for ( j=0; j<M; j++ ) if( j-I ) s += a[i][j] * x[j]; y[i] = ( b[i]-s) / a[i][i]; } max = 0; for ( i=0; i<M; i++ ) if( max < fabs( x[i]-y[i] ) ) max = fabs( x[i] - y[i] ); if ( max < eps ) break; printf( "\nk=%d, ", k ); for( i=0; i<M; i++ ) printf( "y[%d]=%lf ",i, y[i] ); k++; for( i=0; i<M; i++ ) x[i] = y[i]; } if(k == N) { printf( "ERROR!\n" ); return; } printf( "\nk=%d, ", k ); for ( i=0; i < M; i++ ) printf( y[%d]=%lf ", i, y[ i ] ); } 数值分析实验指导书 25 6、〖牛顿法非线性方程求根程序〗 第25页 #include "math.h" double f ( double x ) { return( x - ( x - exp(-x) ) / (1+x) ); } main() { double eps = 1e-10; double x1, x2 = 0.5; do{ x1 = x2; x2 = f( x1 ); printf( "\nx=%20.16lf\n", x2 ); }while( fabs( x2 - x1 ) >= eps ); } 数值分析实验指导书 26 7、〖幂法求矩阵特征值程序〗 第26页 #include "math.h" #define N 3 main() { double a[N][N]={ { - 1, 2, 1},{ 2, - 4, 1 },{ 1, 1, - 6 } }; double s, mu0, mu = -1000, v[N], u[N] = { 1,1,1 }; double eps = 1e-5; int i, j; do{ mu0 = mu; for( i=0; i<N; i++ ) { s = 0; for( j=0; j<N; j++ ) s += a[i][j] * u[j]; v[i] = s; } mu = fabs( v[0] ); for( i=0; i<N; i++ ) if( mu< fabs( v[i] ) ) mu = fabs( v[i] ); for( i=0 ; i<N; i++ ) u[i] = v[i] / mu; }while( fabs( mu-mu0 ) >= eps ); printf( "\nu=%lf\n[", mu ]; for( i=0; i<N; i++ ) printf( "%lf ", u[i] ); printf(") \n" ); } 数值分析实验指导书 27 8、〖四阶经典的龙格-库塔方法解常微分方程程序〗 第27页 float f(float x,float y) { return( 4*x / y – x*y ); } main() { float x0=0,y0=0,x1,y1; float h=0.1,k1,k2,k3,k4; int n=1,N=5; do { x1=x0+h; k1=f(x0,y0); k2=f(x0+h/2,y0+h/2*k1); k3=f(x0+h/2,y0+h/2*k2); k4=f(x1,y0+h*k3); y1=y0+h/6*(k1+2*k2+2*k3+k4); printf("n=%d,x1=%f,y1=%f\n",n,x1,y1); n++; x0=x1;y0=y1; }while(n<N+1); return; }