《数值计算方法实验指导(Matlab版).pdf》由会员分享，可在线阅读，更多相关《数值计算方法实验指导(Matlab版).pdf（38页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、数值计算方法实验指导（Matlab 版）学院 数学与统计学学院计算方法课程组1 / 38数值计算方法实验数值计算方法实验 1 1 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 1算法设计原则验证（之相近数相减、大数吃小数和简化计算步骤）2.2. 实验题目实验题目16 (1) 取z 10，计算z 1 z和1/( z 1 z)，验证两个相近的数相减会造成有效数字的损失 (2) 按不同顺序求一个较大的数（123）与1000 个较小的数（310数吃小数的现象 (3) 分别用直接法和九韶算法计算多项式15）的和，验证大P(x) a0 xn a1xn1 a
2、n1x an在x=1.00037 处的值验证简化计算步骤能减少运算时间对于第(3)题中的多项式P(x)，直接逐项计算需要n (n 1) 21法和n次加法，使用九韶算法n 1次乘2P(x) (a0 x a1)x a2)x an1)x an则只需要n次乘法和n次加法3.3. 实验目的实验目的验证数值算法需遵循的若干规则4.4. 基础理论基础理论设计数值算法时，应避免两个相近的数相减、 防止大数吃小数、简化计算步骤减少运算次数以减少运算时间并降低舍入误差的积累 两相近的数相减会损失有效数字的个数， 用一个大数依次加小数，小数会被大数吃掉，乘法运算次数太多会增加运算时间5.5. 实验环境实验环境操作系
3、统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程 (1) 直接计算并比较； (2) 法 1：大数逐个加 1000 个小数，法 2：先把 1000 个小数相加再与大数加；(3) 将由高次项到低次项的系数保存到数组An中，其中 n 为多项式次数7.7. 结果与分析结果与分析 (1) 计算的z 1分析： (2) 123 逐次加 1000 个310的和是 ，先将 1000 个310相加，再用这个和与 1232 / 3866z= ，1/( z 1 z)相加得分析： (3) 计算次的多项式：直接计算的结果是，用时；用九韶算法计算的结果是，用时分析：8.8. 附录：程序清单附
4、录：程序清单 (1) (1) 两个相近的数相减两个相近的数相减%*%*程 序 名：ex1_1.m *%*程序功能：验证两个相近的数相减会损失有效数字个数 *%*z=1e16;x,y=(2)(2) 大数吃小数大数吃小数%*%*程 序 名：ex1_2.m *%*程序功能：验证大数吃小数的现象. *%*clc; % 清屏clear all; % 释放所有存变量format long; % 按双精度显示浮点数z=123; % 大数t=3e-15; % 小数x=z; % 大数依次加小数3 / 38% 重复 1000 次给 x 中加上 t% x =z 1z；z)；% y =1/( z 1 y=0; % 先
5、累加小数y=z + y; % 再加到大数x,y=(3)(3) 九韶算法九韶算法%*%*程 序 名：ex1_3.m *%*程序功能：验证九韶算法可节省运行时间. *%*clc; % 清屏clear all; % 释放所有存变量format long; % 按双精度显示浮点数A=8,4,-1,-3,6,5,3,2,1,3,2,-1,4,3,1,-2,4,6,8,9,50,-80,12,35,7,-6,42,5,6,23,74,65,55,80,78,77,98,56;A(10001)=0; % 扩展到 10001 项，后面的都是分量 0 % A为多项式系数,从高次项到低次项x=1.00037;n=
6、9000; % n为多项式次数% 直接计算begintime=clock; % 开始执行的时间endtime=clock; % 完毕执行的时间time1=etime(endtime,begintime); % 运行时间disp(直接计算);disp(p(,num2str(x),)=,num2str(p);4 / 38% 累加多项式的 i 次项 % 求 x 的 i 次幂% 重复 1000 次给 y 中加上 tdisp(运行时间: ,num2str(time1),秒);% 九韶算法计算begintime=clock; % 开始执行的时间 % 累加九韶算法中的一项endtime=clock; % 完
7、毕执行的时间time2=etime(endtime,begintime); % 运行时间disp( );disp(九韶算法计算);disp(p(,num2str(x),)=,num2str(p);disp(运行时间: ,num2str(time2),秒);5 / 38数值计算方法实验数值计算方法实验 1 1 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 1算法设计原则验证（之数值稳定性）2.2. 实验题目实验题目计算定积分In其中：x01nex1dx,n 0,1,10，分别用教材例 1-7 推导出的算法 A 和 B，1In1nIn1In1(1
8、In)算法 A：算法 B：nI 0.63210I10 0验证算法不稳定时误差会扩大3.3. 实验目的实验目的验证数值算法需遵循的若干规则4.4. 基础理论基础理论设计数值算法时，应采用数值稳定性好的算法数值稳定的算法， 误差不会放大，甚至会缩小；而数值不稳定的算法会放大误差5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程分别用数组 IA 和 IB 保存两种算法计算的结果7.7. 结果与分析结果与分析运行结果： （或拷屏）n012345678910分析：6 / 38算法 A算法 B精确值8.8. 附录：程序清单附录：程序清单%*%*
9、程 序 名：ex1_4.m *%*程序功能：验证数值稳定性算法可控制误差. *%*clc; % 清屏clear all; % 释放所有存变量format long; % 按双精度显示浮点数I=0.856, 0.144, 0.712, 0.865, .0.538, 0.308, 0.154, 0.938, .0.492, 0.662, 0.843; % 保留 14 位小数的精确值, 是 Matlab 中的续行符% 算法 AIA(1) = 0.6321; % Matlab下标从 1 开始，所以要用 IA(n+1)表示原问题中的 I(n）% 算法 Bdisp(n算法 A算法 B精确值);for n=
10、1:11 fprintf(%2d %14.6f %14.6f %14.6fn,n-1,IA(n),IB(n),I(n);end % n显示为 2 位整数, 其它显示为 14 位其中小数点后显示 6 位的小数7 / 38数值计算方法实验数值计算方法实验 1 1 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 1算法设计原则（除数绝对值不能太小）2.2. 实验题目实验题目将线性方程组增广矩阵利用初等行变换可化为a11a21a12a2221b1r2aa111a11rb20a12b1r1ar a11222a22b20a12b1a22b20由此可解得x1
11、b1/a11,x2 b2/a22分别解增广矩阵为10161 121 2和1621 2101 1的方程组，验证除数绝对值远小于被除数绝对值的除法会导致结果失真3.3. 实验目的实验目的验证数值算法需遵循的若干规则4.4. 基础理论基础理论设计数值算法时， 应避免除数绝对值远小于被除数绝对值的除法， 否则绝对误差会被放大，使结果失真5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程用二维数组 A 和 B 存放方程组的增广矩阵，利用题目所给初等行变换求解方程组7.7. 结果与分析结果与分析第 1 种顺序的方程组的解为x=，y=；第 2 种
12、顺序的方程组的解为x=，y=分析：8.8. 附录：程序清单附录：程序清单%*%*程 序 名：ex1_5.m *%*程序功能：验证除数的绝对值太小可能会放大误差. *%*8 / 38clc;A=1e-16, 1, 1; 2, 1, 2;B=2, 1, 2; 1e-16, 1, 1; % 增广矩阵% 方程组 A % m = - a_21/a_11 是第 2 行加第 1 行的倍数 % 消去 a_21 % m = - a_12/a_22 是第 1 行加第 2 行的倍数 % 消去 a_12, 系数矩阵成对角线 % 未知数 x1 的值 % 未知数 x2 的值disp(方程组 A 的解: x1=,num2s
13、tr(A(1,3), x2=,num2str(A(2,3);disp( );% 方程组 B % m = - b_21/b_11 是第 2 行加第 1 行的倍数 % 消去 b_21 % m = - b_12/b_22 是第 1 行加第 2 行的倍数 % 消去 b_12, 系数矩阵成对角线 % 未知数 x1 的值 % 未知数 x2 的值disp(方程组 B 的解: x1=,num2str(B(1,3), x2=,num2str(B(2,3);9 / 38数值计算方法实验数值计算方法实验 2 2 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 2非线性
14、方程的迭代解法(之简单迭代法)2.2. 实验题目实验题目用简单迭代法求方程x 4x 10 0在区间 1,2的一个实根，取绝对误差限为321043.3. 实验目的实验目的掌握非线性方程的简单迭代法4.4. 基础理论基础理论简单迭代法：将方程f (x) 0改写成等价形式x (x)，从初值x0开始，使用迭代公式xk1(xk)可以得到一个数列，若该数列收敛，则其极限即为原方程的解 取数列中适当的项可作为近似解5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程7.7. 结果与分析结果与分析8.8. 附录：程序清单附录：程序清单数值计算方法实验
15、数值计算方法实验 2 2 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 2非线性方程的迭代解法(之 Newton 迭代法)2.2. 实验题目实验题目用 Newton 迭代法求方程x 4x 10 0在区间1,2的一个实根，取绝对误差限为321043.3. 实验目的实验目的掌握求解非线性方程的 Newton 迭代法4.4. 基础理论基础理论10 / 38Newton 迭代法：解方程f (x) 0的 Newton 迭代公式为xk1 xk5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程7.
16、7. 结果与分析结果与分析8.8. 附录：程序清单附录：程序清单f (xk)f (xk)数值计算方法实验数值计算方法实验 2 2 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 2非线性方程的迭代解法(之对分区间法)2.2. 实验题目实验题目用对分区间法求方程x x1 0在区间1, 1.5的一个实根，取绝对误差限为103.3. 实验目的实验目的掌握求解非线性方程的对分区间法4.4. 基础理论基础理论对分区间法：取a，b的中点p，若f(p) 0或ba ，则p为方程f (x) 0的近似解；若f(a) f(p) 0，则说明根在区间取a，p中；否则，根在
17、区间取p，b中将新的有根区间记为 a1，b1，对该区间不断重复上述步骤，即可得到方程的近似根5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程用宏定义函数f(x)；为了循环方便，得到的新的有根区间始终用a，b表示；由于新的有根区间可能仍以a为左端点，这样会反复使用函数值f(a)，为减少运算次数，将这个函数值保存在一个变量fa 中；同样在判断新的有根区间时用到函数值f(p)，若新的有根区间以p为左端点，则下一次用到的f(a)实际上就是现在的f(p)，为减少运算次数，将这个函数值保存在一个变量fp中算法的伪代码描述：算法的伪代码描述：I
18、nput：区间端点 a，b；精度要求（即误差限）；函数f（x） ；最大对分次数 NOutput：近似解或失败信息行号伪代码11 / 38注释341234567 8 91011121314151617n 1；fa f（a） ；while n N dop （a+b）/ 2；fp f(p)；对分次数计数器左端点的函数值区间中点中点的函数值 iffp=0or(b-a)/2函数值为 0 或半区间长不超 0 thena p；fa fp；elseb p；endifenddoreturn 错误信息输出错误信息并完毕程序输出近似解并退出程序计数器加一若中点与左端点函数值同号新区间取右半区间否则新区间取左右半区间
19、7.7. 结果与分析结果与分析8.8. 附录：程序清单附录：程序清单说明说明: : 源程序中带有数字的空行，对应着算法描述中的行号源程序中带有数字的空行，对应着算法描述中的行号%*%*程 序 名：Bisection.m *%*程序功能：使用二分法求解非线性方程. *%*f=inline(x3-x-1); % 定义函数 f(x)a=input(有根区间左端点: a=);b=input(右端点:b=);epsilon=input(误差限:epsilona=);N=input(最大对分次数: N=); 1 % 对分次数计数器 n 置 1 2 % 左端点的函数值给变量fafprintf(n k p f
20、(p) a(k) f(a(k);fprintf( b(k) b-an);% 显示表头12 / 38fprintf(%2d%36.6f%12.6f%12.6f%12.6fn,0,a,fa,b,b-a); % 占 2 位其中 0 位小数显示步数 0, 共 12 位其中小数 6 位显示各值 3% while n N% 取区间中点 p% 求 p 点函数值给变量 fp% 输出迭代过程中的中点信息p 和 f(p) 4 5 fprintf(%2d%12.6f%12.6f,n,p,fp); 6 % 如果 f(p)=0 或 b-a 的一半小于误差限 fprintf(nn 近似解为:%fn,p);return;%
21、 并完毕程序 (7) 8 9 % 计数器加 1 10% 若 f(a)与 f(p)同号% 则取右半区间为新的求根区间, 即 a 取作 p% 否则 11% 则输出近似根 p (7) 12 % 保存新区间左端点的函数值 13 14 % 左半区间为新的求根区间, 即 b 取作 p 15 fprintf(%12.6f%12.6f%12.6f%12.6fn,a,fa,b,b-a);%显示新区间端点与左端函数值、区间长度 16fprintf(nn 经过%d 次迭代后未达到精度要求.n,N);% 输出错误信息(行 17)数值计算方法实验数值计算方法实验 2 2 报告报告班级：20#级#x 班学号：20#240
22、9#：#X成绩：1.1. 实验名称实验名称实验 2非线性方程的迭代解法(之 Aitken-Steffensen加速法)2.2. 实验题目实验题目用 Aitken-Steffensen加速法求方程x 4x 10 0在区间1,2的一个实根，取绝对误差限为103.3. 实验目的实验目的熟悉求解非线性方程的Aitken-Steffensen 加速法4.4. 基础理论基础理论13 / 38432将方程f (x) 0改写成等价形式x (x)，得到从初值x0开始的迭代公式xk1(xk)后，基于迭代公式xk1(xk)的 Aitken-Steffensen 加速法是通过“迭代-再迭代-加速”完成迭代的，具体过程
23、为yk(xk),zk(yk),xk1(yk zk)2 xkzk2yk xk5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程为了验证 Aitken-Steffensen加速法可以把一些不收敛的迭代加速成迭代收敛， 我们使用将方程组变形为x 1110 x3，取迭代函数(x) 10 x3，并利用宏定义出迭代22函数由于不用保存迭代过程，所以用 x0 表示初值同时也存放前一步迭代的值，y 和 z 是迭代过程中产生的yk和zk，x 存放新迭代的结果算法的伪代码描述：算法的伪代码描述：Input：初值 x0；精度要求（即误差限）；迭代函数（x
24、） ；最大迭代次数 NOutput：近似解或失败信息行号123456 7 8 9101112n 1；while n N doy （x0） ；z （y） ；2伪代码迭代再迭代注释迭代次数计数器x x0 (y-x0) / (z - 2 y + x0)加速if | x x0| thenreturn x； endifn n + 1；x0 x；enddoreturn 错误信息输出错误信息并完毕程序如果达到精度要求则输出近似值并退出程序计数器加一新近似值给 x0做下次的初值7.7. 结果与分析结果与分析8.8. 附录：程序清单附录：程序清单%*%*程 序 名：Aitken_Steffensen.m *%*
25、程序功能：用 Aitken-Steffensen加速法求方程. *%*clc;clear all;14 / 38phi=inline(0.5 * sqrt( 10 - x3); % 迭代函数x0=input(初值: x0 = );epsilon=input(误差限: epsilon=);N=input(最大迭代次数: N=);disp( n迭代中间值 y(n-1)再迭代结构 z(n-1)加速后的近似值 x(n);fprintf(%2d%54.6fn,0,x0); % 占 2 位整数显示步数 0, 为了对齐, 占 54 位小数 6 位显示 x0 1 % n是计数器 2 % while n K t
26、hen return False；endifn 1；F0 f(x0)；while n N doF0 f(x0)；if F0= 0 then return False；endiflambda 1；k 1；while k K dox x0lambda*F0 / F0Fx f(x)；if |Fx| |F0| then伪代码若该点导数为 0则无法进行迭代，完毕程序下山因子从 1 试起尝试次数计数器Newton 下山公式判断下山条件退出尝试循环下山因子减半尝试次数计数器加 1如果因超过尝试次数退出循环则提示错误信息并完毕程序注释迭代次数计数器16 / 382223242526272829if | x x
27、0| thenreturn x； endif x0 x；F0 Fx；n n + 1；enddoreturn 错误信息否则时尝试成功退出上边循环如果达到精度要求则输出近似值并退出程序新近似值给 x0做下次的初值所求函数值下次也用到计数器加一输出错误信息并完毕程序7.7. 结果与分析结果与分析8.8. 附录：程序清单附录：程序清单%*%*程 序 名：NewtonDownhill.m *%*程序功能：用 Newton 下山法求解非线性方程. *%*clc;clear all;f=inline(x3-x-1); % 函数 f(x)df=inline(3*x2-1); % 函数 f(x)的导函数x0=i
28、nput(初值: x0 = );epsilon=input(误差限: epsilon=);N=input(最大迭代次数: N=);K=input(最大下山尝试次数: K=); 1 % 迭代次数计数器 2 % 存 x0 点函数值fprintf(nn n x(n) f(x(n)n);% 显示表头fprintf(%2d%14.6f%14.6fn,0,x0,fx0);% 2 位整数显示0, 共 14 位小数 6 位显示x0 和 fx0 3 % while n N disp();% 换行显示下山尝试过程的表头 disp(下山因子尝试 x(n)对应 f(x(n)满足下山条件);disp(); 4% 存 x
29、0 点导数值, 每次下山尝试不用重新计算17 / 38ifdfx0=0% 导数为 0 不能迭代disp(无法进行 Newton 迭代);return; end lambda=1.0; % 下山因子从 1 开始尝试k=1; % k下山尝试次数计数器while k=K % 下山最多尝试 K 次% 下山公式 fx=f(x);% 函数值 fprintf(%22.6f%14.6f%14.6f,lambda,x,fx);% 显示尝试结果if (abs(fx)K fprintf(n下山条件无法满足, 迭代失败.nn); return; endfprintf(%2d%14.6f%14.6fn,n,x,fx);
30、% 2 位整数显示步数 n, 共 14 位小数 6 位显示下步迭代结果 22% 达到精度要求否 fprintf(nn方程的近似解为: x%fnn,x);% (23) return;% 达到, 则显示结果并完毕程序(23) end% (24) % 用 x0,fx0 存放前一步的近似值和它的函数值, 进行循环迭代 25 26 27 28fprintf(n迭代%d 次还不满足误差要求.nn,N);18 / 38数值计算方法实验数值计算方法实验 2 2 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 2非线性方程的迭代解法（之弦截法）2.2. 实验题目实
31、验题目用弦截法求方程x 4x 10 0在区间1,2的一个实根，取绝对误差限为103.3. 实验目的实验目的熟悉非线性方程的弦截法4.4. 基础理论基础理论将 Newton 迭代法中的导数用差商代替，得到弦截法(或叫正割法)公式324xk1 xkxk xk1f (xk)f (xk) f (xk1)5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程不保存迭代过程， 所以始终以x0和x1分别存放xk-1和xk， 而 x 存放新产生的迭代值xk+1，这样，下一次迭代时需要把上一步的x1（即xk）赋值于x0（做新的xk-1） 这些点的函数值会
32、重复用到，在迭代公式中也要用到，上一步的x1作为下一步的x0也会再一次用它的函数值，为减少重新计算该点函数值的运行时间，将x1点的函数值保存在变量 fx1 中算法的伪代码描述：算法的伪代码描述：Input：初值 x0，x1；精度要求（即误差限）；函数f（x） ；最大迭代次数 NOutput：近似解或失败信息行号123456 7 8 9101112伪代码fx0 f(x0)；fx1 f(x1)；n 2；while n N do迭代次数计数器(从 2 开始)注释初值点的函数值x x1(x1-x0)*f(x1) / (f(x1) 迭代f(x0)；if | x x0| j % 在下三角部分, 则取 A
33、对于的元素显示 fprintf( %8.4f,A(i,j); elseif i=j % 在对角线上, 则显示 1 fprintf( %8d,1); else % 在上三角部分, 则显示 0 fprintf( %8d,0); end end fprintf(n); % 换行enddisp(U=);for i=1:n for j=1:n if i=j % 在下三角部分或对角线上, 则取 A 对于的元素显示fprintf(%10.4f,A(i,j);% 计算其他的l_ik% 和l_ki else % 在上三角部分, 则显示 0 fprintf(%10d,0); end end fprintf(n);
34、 % 换行end数值计算方法实验数值计算方法实验 3 3 报告报告班级：20#级#x 班：1.1. 实验名称实验名称实验 3解线性方程组的直接法(之改进的 Cholesky 分解)2.2. 实验题目实验题目对矩阵 A 进行改进的 Cholesky 分解，其中#X学号：20#2409#成绩： 62A 113.3. 实验目的实验目的理解矩阵改进的 Cholesky 分解4.4. 基础理论基础理论141014101321矩阵的改进的Cholesky分解是指将矩阵A (aij)nn可以分解为一个单位下三角矩阵L和对角矩阵D与L转置的乘积，即A=LDL，其中L和D各元素可依如下顺序公式计算T30 / 3
35、8k12d ad lkkkr krr1k1lik (aikdtlitlkt)/dk,t1i k 1,k 2,n其中k = 1，2，n5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：VC+6.6. 实验过程实验过程（1）按计算公式依次先计算D的一个元素dk，再计算L中这列的元素lik，且对称位置的元素也取同一个值；（2）因为计算完dk和lij后，akk或aij就不再使用，为节省存储空间，将计算的akk或lij仍存放在矩阵 A 中的相应位置；（3）使用 L 矩阵时需要根据元素所在位置取固定值或A 中相应位置的值L 对角线和上三角部分为 0，下三角部分为A 中对应的元素；D
36、 对角线为 A 中对应的元素，其余都是0算法的伪代码描述：算法的伪代码描述：Input：阶数 n；矩阵 AOutput：矩阵 L (合并存储在数组 A 中)行号12for k 1 to n2dk akkdrlkrr1k1伪代码注释3456 77.7. 结果与分析结果与分析for i k to nk1t1lik (aikdtlitlkt)/dkendforendforreturn L, D计算结果存放在aij输出 L 和 D8.8. 附录：程序清单附录：程序清单%*%*程 序 名：ImprovedCholesky.m *%*程序功能：对称正定矩阵的改进的Cholesky 分解. *%*31 /
37、38n=4; % 矩阵阶数A=6,2,1,-1; 2,4,1,0; 1,1,4,-1; -1,0,-1,3;disp(A =);for i=1:nfor j=1:nfprintf(%10.4f,A(i,j); % 共占 14 位% 一行完毕换行 end fprintf(n);end% Cholesky 分解for k=1:n% 计算 D 对角线上的u_kk% 计算 L 的元素l_ik% 和 L 转置的元素l_kiend % L在 A 下三角, D 在对角线disp(分解结果：);disp(L=);for i=1:nfor j=1:nif ij % 在下三角部分, 则取 A 对于的元素显示fpr
38、intf(%10.4f,A(i,j); elseif i=j % 在对角线上, 则显示 1 fprintf(%10d,1); else % 在上三角部分, 则显示 0 fprintf(%10d,0); end end fprintf(n); % 换行end32 / 38disp(D=);for i=1:nfor j=1:nif i=j % 在对角线上, 则取 A 对于的元素显示fprintf(%10.4f,A(i,j);else % 其余显示 0 fprintf(%10d,0); end end fprintf(n); % 换行end数值计算方法实验数值计算方法实验 3 3 报告报告班级：20
39、#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 3解线性方程组的直接法(之追赶法)2.2. 实验题目实验题目用追赶法求解线性方程组2100 1320 x10 x120121x03 0035x413.3. 实验目的实验目的熟悉解线性方程组的追赶法4.4. 基础理论基础理论对于系数矩阵为三对角矩阵的方程组，其Crout 分解可分解为b1c1a s12b2c2a2s1t11tA 22can1bn1n1an1sn11anbnansn这样，解方程组可以由如下2 步完成：“追”：s1 b1,y1 f1/s1,ti1 ci1/si1,si biaiti1,yi ( fiaiyi
40、1)/si,i 2,3,n,33 / 38tn11T其中：( f1, fn)为方程组的常数项，tn没用；“赶”：xn yn,xi yitixi1,i n1,n2,1.5.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程在“追”的过程中，向量s和 y 都有n个元素，t只有n-1 个元素，又s1和y1的计算公式与其它si和yi不同，所以先单独计算s1和y1，然后在一个n-1 次循环中，求其它si和yi以与ti由于在“追”的过程中，bi，ci和fi在分别计算完对应的si，ti和yi后就不再使用，所以借用数组b，c和f存储向量s，t和y；同样
41、在“赶”的过程中，yi在计算完对应的xi后就不再使用，所以再一次借用数组f 存储向量x追赶法算法的伪代码描述：追赶法算法的伪代码描述：Input：阶数 n；三对角矩阵的三条对角线向量a，b，c，常数项向量 fOutput：方程组的解 x行号123456 7 8 9 10 11 12改进的追赶法算法的伪代码描述：改进的追赶法算法的伪代码描述：Input：阶数 n；三对角矩阵的三条对角线向量a，b，c，常数项向量 fOutput：方程组的解 x行号1s1 b1；34 / 38伪代码注释s 向量仍存放 b 中s1 b1；y1 f1/s1；for k 2 to ntk-1 ck-1/sk-1；sk b
42、k - ak tk-1； yk (fk - ak yk-1) / sk；endforxn yn；for k n-1 to 1 xk yk - tk yk+1；endforreturn x赶伪代码追注释23456 7 8 9 10 11 12f1 f1/b1；for k 2 to nck-1 ck-1/sk-1；bk bk - ak ck-1； fk (fk - ak fk-1) / bk；endforxn fn；for k n-1 to 1 fk fk - ck fk+1；endforreturn fy 向量暂存在 f 中t 向量仍存放才中x 向量也存在 f 中7.7. 结果与分析结果与分析8
43、.8. 附录：程序清单附录：程序清单%*%*程 序 名：ChaseAfter.m *%*程序功能：用追赶法求解三对角线性方程组. *%*clc;clear all;n=4;a=0,-1,-1,-3;b=2, 3, 2, 5;c=-1, -2, -1, 0;f=0, 1, 0, 1;% 追s(1) = b(1);y(1) = f(1); % 先单独求 s_1 和 y_1for k = 1 : n-1% 再求 t_i(i=1,2,n-1)% s_i(i=2,3,n)% y_i(i=2,3,n)end% 赶x(n) = y(n); % 先单独求 x_nfor k = n-1 : -1 : 1% 再求
44、 x_i(i=n-1,n-2,1)35 / 38endx=x % 输出解向量-改进的程序:%*%*程 序 名：ChaseAfter.m *%*程序功能：用追赶法求解三对角线性方程组. *%*clc;clear all;n=4;a=0,-1,-1,-3;b=2, 3, 2, 5;c=-1, -2, -1, 0;f=0, 1, 0, 1;% 追% b(1)=b(1);% s_1 仍在 b_1 中，不用重新计算y(1)=f(1)/b(1);% 先单独 y_1for k=1:n-1% 再求 t_i(i=1,2,n-1)% s_i(i=2,3,n)% y_i(i=2,3,n)end% 赶% f(n)=f
45、(n);for k=n-1:-1:1endx=f% x_n 等于 y_n 仍在 f_n 中% 再求 x_i(i=n-1,n-2,1)% 输出解向量数值计算方法实验数值计算方法实验 4 4 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 4解线性方程组的迭代法(之 Jacobi 迭代)2.2. 实验题目实验题目用 Jacobi 迭代法求解线性方程组36 / 38 x1x12x12x2x22x22x3x3x3513任取3.3. 实验目的实验目的掌握解线性方程组的 Jacobi 迭代法4.4. 基础理论基础理论将第i（1 i n）个方程ai1x1 ai
46、2x2 ainxn bi移项后得到等价方程xi (bi ai1x1 ai,i1xi1 ai,i1xi1 ainxn)/aii便可构造出 Jacobi 迭代公式k)(k)(k)xi(k1) (biai1x1(k)ai,i1xi(,1ai,i1xi1ainxn)/aii,k 0,15.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程7.7. 结果与分析结果与分析8.8. 附录：程序清单附录：程序清单数值计算方法实验数值计算方法实验 4 4 报告报告班级：20#级#x 班学号：20#2409#：#X成绩：1.1. 实验名称实验名称实验 4解
47、线性方程组的迭代法(之 Gauss-Seidel 迭代)2.2. 实验题目实验题目用 Gauss-Seidel 迭代法求解线性方程组37 / 3810 x1 x12x1 x211x2 x23x2 2x3 x310 x3 x33x4 x48x462511153.3. 实验目的实验目的掌握解线性方程组的 Gauss-Seidel 迭代法4.4. 基础理论基础理论将第i（1 i n）个方程ai1x1 ai2x2 ainxn bi移项后得到等价方程xi (bi ai1x1 ai,i1xi1 ai,i1xi1 ainxn)/aii便可构造出 Jacobi 迭代公式(k)k)(k)(k)xi(k1) (biai1x1ai,i1xi(,1 i,i1xi1ainxn)/aii,k 0,1和 Gauss-Seidel 迭代公式k1)k)(k)xi(k1) (biai1x1(k1)i,i1xi(ai,i1xi(,11ainxn)/aii,k 0,15.5. 实验环境实验环境操作系统：Windows xp； 程序设计语言：Matlab6.6. 实验过程实验过程7.7. 结果与分析结果与分析8.8. 附录：程序清单附录：程序清单38 / 38