数值计算方法实验指导讲解

《数值计算方法》

实验指导（Matlab版）肇庆学院数学与统计学学院

计算方法课程组《数值计算方法》实验1报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验1算法设计原则验证(之相近数相减、大数吃小数和简化计算步骤)实验题目取z=1。16,计算-云和1/(t7军1+云)，验证两个相近的数相减会造成有效数字的损失.按不同顺序求一个较大的数(123)与1000个较小的数(3X10-15)的和，验证大数吃小数的现象.分别用直接法和秦九韶算法计算多项式P(x)=axn+axn-1++ax+a在x=1.00037处的值.验证简化计算步骤能减少运算时间.n+1，一,对于第(3)题中的多项式P(x),直接逐项计算需要n+(n-1)+•••+2+1=——次乘法和n次加法，使用秦九韶算法P(x)=(((ax+a)x+a)x—+a)x+a则只需要n次乘法和n次加法.实验目的验证数值算法需遵循的若干规则.基础理论设计数值算法时，应避免两个相近的数相减、防止大数吃小数、简化计算步骤减少运算次数以减少运算时间并降低舍入误差的积累.两相近的数相减会损失有效数字的个数，用一个大数依次加小数，小数会被大数吃掉，乘法运算次数太多会增加运算时间.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程直接计算并比较；法1：大数逐个加1000个小数，法2：先把1000个小数相加再与大数加；将由高次项到低次项的系数保存到数组A[n]中，其中n为多项式次数.结果与分析(1)计算的k：z+1—z=，1/(yZ+1+tz).分析：(2)123逐次加1000个3x10-6的和是，先将1000个3X10-6相加，再用这个和与123相加得.分析：(3)计算次的多项式：直接计算的结果是，用时;用秦九韶算法计算的结果是，用时分析：8.附录：程序清单两个相近的数相减.%**个**********************************************************%*程序名：ex1_1.m*%*程序功能：验证两个相近的数相减会损失有效数字个数%*************************************************************z=1e16;%x=气Z+1—tZ；%y=1/"E+.&)x,y(2)大数吃小数%*************************************************************%*程序名：ex1_2.m*%*程序功能：验证大数吃小数的现象.%*************************************************************clc;%清屏clearall;%释放所有内存变量formatlong;%按双精度显示浮点数z=123;%大数t=3e-15;%小数x=z;%大数依次加小数%重复1000次给x中加上t

y=0；%先累加小数%重复1000y=0；%先累加小数%重复1000次给y中加上ty=z+y;%再加到大数(3)秦九韶算法%**个**********************************************************%*程序名：ex1_3.m*%*程序功能：验证秦九韶算法可节省运行时间.%*************************************************************clc;%清屏clearall;%释放所有内存变量formatlong;%按双精度显示浮点数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=9000;%n为多项式次数%直接计算begintime=clock;P=0:鼬i=「：-l』T=1:奴k=l:it=t*京.endp=p-A(n-i-L)%;end.endtime=clock;time1=etime(endtime,begintime);disp('直接计算');disp(['p(',num2str(x),')=',num2str(p)]);%开始执行的时间%求x的i次幂%累加多项式的i次项%结束执行的时间%运行时间disp(['运行时间：',num2str(time1),'秒秒]);

%%开始执行的时间%求x的i次幂%累加多项式的i次项%结束执行的时间%运行时间%秦九韶算法计算%开始执行的时间%累加秦九韶算法中的一项begintime=clock;fori=0;p-=p-A{i-1):endendtime=clock;disp(['运行时间：',num2str(time2),秒]);《数值计算方法》实验1报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验1算法设计原则验证（之数值稳定性）实验题目计算定积分I=j1心心一i&,n=0,1,10，分别用教材例1-7推导出的算法A和B，n0其中：算法A：算法A：<nI0=1一nIn-1r0.63216II一(1-I)算法B：'n-1nn110^0验证算法不稳定时误差会扩大.实验目的验证数值算法需遵循的若干规则.基础理论设计数值算法时，应采用数值稳定性好的算法.数值稳定的算法，误差不会放大，甚至会缩小；而数值不稳定的算法会放大误差.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程分别用数组IA[]和IB[]保存两种算法计算的结果.8.附录：程序清单%**个**********************************************************%*程序名：ex1_4.m%*程序功能：验证数值稳定性算法可控制误差.%*************************************************************clc;%清屏clearall;%释放所有内存变量formatlong;%按双精度显示浮点数I=[0.63212055882856,0.36787944117144,0.26424111765712,0.20727664702865,0.17089341188538,0.14553294057308,0.12680235656154,0.11238350406938,.0.10093196744492,0.09161229300662,0.08387707010843];%保留14位小数的精确值，•是Matlab中的续行符%算法AIA(1)=0.6321;%Matlab下标从1开始，所以要用IA(n+1)表示原问题中的I(n)forn=1:10IA(n-L)=1-n*lA(n):end.%算法B0:forn=10:-1:1/n:enddisp('n算法A算法B精确值')；forn=1:11fprintf('%2d%14.6f%14.6f%14.6f\n',n-1,IA(n),IB(n),I(n));end%n显示为2位整数，其它显示为14位其中小数点后显示6位的小数《数值计算方法》实验1报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验1算法设计原则（除数绝对值不能太小）实验题目将线性方程组增广矩阵利用初等行变换可化为(a11a12a22r—%r2an1T01a12a22r—―r1。22Tf(a11a12a22r—%r2an1T01a12a22r—―r1。22Tfa11〔0由此可解得T叩an，%分别解增广矩阵为实验目的验证数值算法需遵循的若干规则.基础理论设计数值算法时，应避免除数绝对值远小于被除数绝对值的除法，否则绝对误差会被放大，使结果失真.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程用二维数组A和B存放方程组的增广矩阵，利用题目所给初等行变换求解方程组.结果与分析第1种顺序的方程组的解为x=，y=；第2种顺序的方程组的解为x=，y=.分析：8.附录：程序清单%**个**********************************************************%*程序名：ex1_5.m*%*程序功能：验证除数的绝对值太小可能会放大误差.%*************************************************************

clc;A=[1e-16,1,1;2,1,2];B=[2,1,2;1e-16,1,1];%方程组A1);m=-A(l:2)/A(2;2);A(1:：)=A(1:A(l;3)=A(L3)/A(l=1);%增广矩阵%m=-a_{21}/a_{11}是第2行加第1行的倍数%消去a_{21}%m=-a_{12}/a_{22}是第1行加第2行的倍数%消去a_{12},系数矩阵成对角线%未知数x1的值%未知数x2的值disp(['方程组A的解:x1=',num2str(A(1,3)),',x2=',num2str(A(2,3))]);disp('');%方程组B%m=-b_{21}/b_{11}是第2行加第1clc;A=[1e-16,1,1;2,1,2];B=[2,1,2;1e-16,1,1];%方程组A1);m=-A(l:2)/A(2;2);A(1:：)=A(1:A(l;3)=A(L3)/A(l=1);%增广矩阵《数值计算方法》实验2报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验2非线性方程的迭代解法(之简单迭代法)实验题目用简单迭代法求方程尤3+4x2-10=0在区间［1,2］内的一个实根，取绝对误差限为10-4.实验目的掌握非线性方程的简单迭代法.基础理论简单迭代法：将方程f(x)=0改写成等价形式x=9(x)，从初值x0开始，使用迭代公式x*=9(xk)可以得到一个数列，若该数列收敛，则其极限即为原方程的解.取数列中适当的项可作为近似解.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程结果与分析附录：程序清单《数值计算方法》实验2报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验2非线性方程的迭代解法(之Newton迭代法)实验题目用Newton迭代法求方程尤3+4尤2-10=0在区间［1,2］内的一个实根，取绝对误差限为10-4.实验目的掌握求解非线性方程的Newton迭代法.基础理论Newton迭代法：解方程/(x)=0的Newton迭代公式为x=x-k.Ef'(x)k实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程结果与分析附录：程序清单《数值计算方法》实验2报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验2非线性方程的迭代解法(之对分区间法)实验题目用对分区间法求方程x3-x-1=0在区间［1,1.5］内的一个实根，取绝对误差限为10-4.实验目的掌握求解非线性方程的对分区间法.基础理论对分区间法：取［sb］的中点p,若f(p)=0或b-a<们则p为方程f(x)=0的近似解；若f(a)f(p)<0，则说明根在区间取［a,p］中；否则，根在区间取［p,b］中.将新的有根区间记为［a”b1］,对该区间不断重复上述步骤，即可得到方程的近似根.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程用宏定义函数f(x)；为了循环方便，得到的新的有根区间始终用［a，b］表示；由于新的有根区间可能仍以a为左端点，这样会反复使用函数值f(a)，为减少运算次数，将这个函数值保存在一个变量fa中；同样在判断新的有根区间时用到函数值f(p)，若新的有根区间以p为左端点，则下一次用到的f(a)实际上就是现在的f(p)为减少运算次数，将这个函数值保存在一个变量fp中.算法的伪代码描述：Input:区间端点a，b；精度要求(即误差限)&函数f(x)；最大对分次数NOutput:近似解或失败信息行号伪代码注释1n—1；对分次数计数器2fa—f(a)；左端点的函数值3whilen<Ndo4p—(a+b)/2；区间中点5fp—Hp)；中点的函数值6iffp=0or(b-a)/2<£then函数值为0或半区间长不超£7returnp；输出近似解并退出程序8endif9n—n+1；计数器加一10iffa-fp>0then若中点与左端点函数值同号11a—p；新区间取右半区间12fa—fp；13else否则14b—p；新区间取左右半区间

15endif16enddo17return错误信息输出错误信息并结束程序结果与分析附录：程序清单说明：源程序中带有数字的空行，对应着算法描述中的行号rvf““““““““““““““““““““““““““““““““““““““““““““““““““““““““““%**个*******************************************************%*程序名：Bisection.m*%*程序功能：使用二分法求解非线性方程.rvf““““““““““““““““““““““““““““““““““““““““““““““““““““““““““%**********************************************************f=inline('xAf=inline('xA3-x-1');%定义函数f(x)a=input('有根区间左端点：a=');b=input('右端点:b=');epsilon=input('误差限:epsilona=');N=input('最大对分次数:N=');%对分次数计数器n置%对分次数计数器n置1%左端点的函数值给变量fafprintf('\nkpfprintf('b(k)f(p)a(k)f(a(k))');b-a\n');%显示表头fprintf('%2d%36.6f%12.6f%12.6f%12.6f\n',0,a,fa,b,b-a);%占2位其中0位小数显示步数0,共12位其中小数6位显示各值3%whilen<N4%取区间中点p5%求p点函数值给变量fpfprintf('%2d%12.6f%12.6f,n,p,fp);%输出迭代过程中的中点信息p和f(p)6fprintf('\n\n近似解为:％f\n',p);return;8%如果f(p)-0或b-a的一半小于误差限£%则输出近似根p(7)%并结束程序(7)%计数器加1%若f(a)与f(p)同号11%则取右半区间为新的求根区间，即a取作p12%保存新区间左端点的函数值13%否则14%左半区间为新的求根区间，即b取作p15fprintf('%12.6f%12.6f%12.6f%12.6f\n',a,fa,b,b-a);%显示新区间端点及左端函数值、区间长度16fprintf('\n\n经过％d次迭代后未达到精度要求.\n',N);%输出错误信息(行17)《数值计算方法》实验2报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验2非线性方程的迭代解法(之Aitken-Steffensen加速法)实验题目用Aitken-Steffensen加速法求方程尤3+4尤2—10=0在区间［1,2］内的一个实根，取绝对误差限为10-4.实验目的熟悉求解非线性方程的Aitken-Steffensen加速法.基础理论将方程f(x)=0改写成等价形式x=9(x)，得到从初值七开始的迭代公式气+］=甲(xk)后，基于迭代公式xk+1=9(xk)的Aitken-Steffensen加速法是通过“迭代-再迭代-加速”完成迭代的，具体过程为J=9(x),乙=93),x=x—(七七*.

kkkkk+1kz—2y+x实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程为了验证Aitken-Steffensen加速法可以把一些不收敛的迭代加速成迭代收敛，我们使用将方程组变形为x=1<10—x3，取迭代函数9(x)=；(10—x3，并利用宏定义出迭代函数.由于不用保存迭代过程，所以用x0表示初值同时也存放前一步迭代的值，y和z是迭代过程中产生的耻和zk，x存放新迭代的结果.算法的伪代码描述：Input：初值x0；精度要求(即误差限)&迭代函数p(x)；最大迭代次数NOutput:近似解或失败信息行号伪代码注释1n—1；迭代次数计数器2whilen<Ndo3y—p(x0)；迭代4z—p(y)；再迭代5x—x0—(y-x0)2/(z-2y+x0)加速6if1x一x0|<£then如果达到精度要求7returnx；则输出近似值并退出程序8endif

9n—n+1；计数器加一10x0—x；新近似值给x0做下次的初值11enddo12return错误信息输出错误信息并结束程序7.结果与分析8.附录：程序清单rvf“““““““““““““““““““““““““““““““““““““““““““““““““““““““““““““%**个**********************************************************%*程序名：Aitken_Steffensen.m%*程序功能：用Aitken-Steffensen加速法求方程.rvf“““““““““““““““““““““““““““““““““““““““““““““““““““““““““““““%*************************************************************clc;clearall;phi=inline('0.5*sqrt(10-xA3)');%迭代函数x0=input(初值:x0=');epsilon=input('误差限：epsilon=');N=input('最大迭代次数:N=');disp('n迭代中间值y(n-1)再迭代结构z(n-1)加速后的近似值x(n)');fprintf('%2d%54.6f\n',0,x0);%占2位整数显示步数0,为了对齐，占54位小数6位显示x0%n是计数器%n是计数器%whilen<=N%迭代%再迭代%加速z是中间变量,x是下一步的近似值y=3;z=3;x=3;%x0初值及前一步的近似值,y和fprintf('%2d%18.6f%18.6f%18.6f\n',n,y,z,x);%显示中间值和迭代近似值6%如果与上一步近似解差的绝对值不超过误差限fprintf('\n\n近似解x~x(%d)~%f\n',n,x);%则输出近似根(7)，可简略为:fprintf('\n\n近似解x=%f,x);return;%并结束程序(7)8%相当于endif9%计数器加110%新近似值x作为下一次迭代的初值11fprintf('\n迭代％d次还不满足误差要求.\n\n',N);%输出错误信息(12)《数值计算方法》实验2报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验2非线性方程的迭代解法(之Newton下山法)实验题目用Newton下山法求方程尤3+4尤2-10=0在区间［1,2］内的一个实根，取绝对误差限为10-4.实验目的熟悉非线性方程的Newton下山法.基础理论Newton下山法：Newton下山法公式为^=^-人H「，使I/(x)1<1/(x)1,k+1kk广(x)k+1kk其中0<Xk<1.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程定义函数f(x)和df(x)，其中df(x)是f(x)的导函数.每步迭代时先取下山因子为1,尝试迭代，判断尝试结果是否满足下山因子，若满足则作为这步的迭代结果；否则将下山因子减半，然后再尝试.为防止当前的xk是极小值点，附近不会有满足下述条件的其它点，使尝试陷入死循环，同时计算机中能表示出的浮点数也有下界，因此我们设置了最大尝试次数.当超过最大尝试次数时，不再进行下山尝试.由于反复尝试迭代且要判断下山条件，所以f(x0)和f(x0)会反复使用，为避免重复计算浪费运行时间，将这两个值分别保存在变量fx0和dfx0.而尝试产生的节点，判断下山条件时要用到它的函数值，若尝试成功，这个点会作为下一步的初值再使用，所以把该点的函数值也保存在变量fx中.算法的伪代码描述：Input:初值x0；精度要求(即误差限)&函数及其导函数f(x)和f(x)；最大迭代次数N；K下山尝试最大次数Output:近似解或失败信息行号伪代码注释1n—1；迭代次数计数器2F0—f(x0)；3whilen<Ndo4F0’—g)；5ifF0’=0then若该点导数为06returnFalse；则无法进行迭代，结束程序7endif

8lambda—1；下山因子入从1试起9k—1；尝试次数计数器10whilek<Kdo11x—x0-lambda*F0/F0’Newton下山公式12Fx—f（x）；13ifIFxl<IF0Ithen判断下山条件14break；退出尝试循环15endif16lambda—lambda/2；下山因子减半17k—k+1；尝试次数计数器加118endwhile19ifk>Kthen如果因超过尝试次数退出循环20returnFalse；则提示错误信息并结束程序21endif否则时尝试成功退出上边循环22ifIx一x0|<£then如果达到精度要求23returnx；则输出近似值并退出程序24endif25x0—X；新近似值给x0做下次的初值26F0—Fx；所求函数值下次也用到27n—n+1；计数器加一28enddo29return错误信息输出错误信息并结束程序7.结果与分析8.附录：程序清单rvf“““““““““““““““““““““““““““““““““““““““““““““““““““““““““““““%**个**********************************************************%*程序名：NewtonDownhill.m*%*程序功能：用Newton下山法求解非线性方程.m“““““““““““““““““““““““““““““““““““““““““““““““““““““““““““““%*************************************************************clc;clearall;%函数f(x)%函数f(x)的导函数f=inline（'xA3-x-1'）;df=inline（'3*xA2-1'）;x0=input（初值:x0='）;epsilon=input（'误差限：epsilon='）;N=input（'最大迭代次数%函数f(x)%函数f(x)的导函数迭代次数计数器存x0点函数值fprintf('\n\nnx(n)f(x(n))\n');fprintf('\n\nnx(n)f(x(n))\n');显示表头fprintf('%2d%14.6f%14.6f\n',0,x0,fx0);%2位整数显示0,共14位小数6位显示x0和fx03%whilen<Ndisp('');%换行显示下山尝试过程的表头disp('下山因子尝试x(n)对应f(x(n))满足下山条件');disp('');4%存x0点导数值，每次下山尝试不用重新计算ifdfx0==0%导数为0不能迭代disp(‘无法进行Newton迭代’)；return;endlambda=1.0;%下山因子从1开始尝试k=1;%k下山尝试次数计数器whilek<=K%下山最多尝试K次k=xO-lambda*fftO'dhtO:%下山公式fx=f(x);fprintf('%22.6f%14.6f%14.6f,lambda,x,fx);%函数值%显示尝试结果if(abs(fx)<abs(fx0))fprintf('满足4');%判断是否满足下山条件break;elsefprintf('不满足4');end%是，则退出下山尝试的循环lambda=lambda/2;%不是，则下山因子减半k=k+1;endifk>K%计数器加1fprintf('\n下山条件无法满足，迭代失败.\n\n');return;endfprintf('%2d%14.6f%14.6f\n',n,x,fx);%2位整数显示步数n,共14位小数6位显示下步迭代结果22%达到精度要求否fprintf('\n\n方程的近似解为：x^f\n',x);%(23)return;%达到，则显示结果并结束程序(23)22end%(24)%用x0,fx0存放前一步的近似值和它的函数值，进行循环迭代25262728fprintf('\n迭代％d次还不满足误差要求.\n\n',N);《数值计算方法》实验2报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验2非线性方程的迭代解法（之弦截法）实验题目用弦截法求方程X3+4X2-10=0在区间［1,2］内的一个实根，取绝对误差限为10-4.实验目的熟悉非线性方程的弦截法.基础理论将Newton迭代法中的导数用差商代替，得到弦截法（或叫正割法）公式X=X一-k__-^f(X).*kf(X)-f(X)kkk-1实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程不保存迭代过程，所以始终以X0和X1分别存放Xk-1和Xk，而x存放新产生的迭代值Xk+1，这样，下一次迭代时需要把上一步的X1（即Xk）赋值于X0（做新的Xk-1）.这些点的函数值会重复用到，在迭代公式中也要用到，上一步的X］作为下一步的X0也会再一次用它的函数值，为减少重新计算该点函数值的运行时间，将X］点的函数值保存在变量fx1中.算法的伪代码描述：Input:初值x0，x1；精度要求（即误差限）£；函数f（x）；最大迭代次数NOutput:近似解或失败信息行号伪代码注释1fx0—f(x0)；初值点的函数值2fx1—f(x,)；3n—2；迭代次数计数器（从2开始）4whilen<Ndo5x—x1-(x1-x0)*f(x1)/(f(x,)-f(x0))；迭代6if1x一x0|<£then如果达到精度要求7returnx；则输出近似值并退出程序8endif9x0—x1；（按顺序）最新的2个近似值，10x1—x；分别做下次的初值x0，x111fx0—fx1；下次初值点的函数值12fx1—f(x1)；13n—n+1；14enddo15return错误信息输出错误信息并结束程序7.结果与分析8.附录：程序清单%**个**********************************************************%*程序名：SecantMethod.m*%*程序功能：用弦截法求解非线性方程.%*************************************************************clc;clearall;f=inline('2*xA3-5*x-1');%函数f(x)x0=input('第一初值:x0=');x1=input('第二初值:x1=');epsilon=input('误差限：epsilon=');N=input('最大迭代次数:N=');fprintf('\nnx(n)\n');fprintf('%2d%14.6f\n',0,x0);fprintf('%2d%14.6f\n',1,x1);%显示表头%占2位显示步数0,共14位其中小数6位显示x0%占2位显示步数1,共14位其中小数6位显示x11%存x0点函数值2%存x1点函数值3%迭代计数器4%whilen<N:―"三土顼%弦截法公式fprintf（'%2d%14.6f\n',n,x）;%显示迭代过程6%达到精度要求否fprintf（'\n\n方程的近似解为:x~%f\n\n',x）;return;%达到，则显示结果并结束程序89%原x1做x0为前两步的近似值W%现x做x1为一两步的近似值11%x0点函数值12%计算x1点函数值，为下一次循环13%计数器加114fprintf（'\n迭代％d次还不满足误差要求.\n\n',N）;

《数值计算方法》实验3报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验3解线性方程组的直接法（之Gauss消去法）实验题目用Gauss消去法求解线性方程组'0.0012.000-1.000'0.0012.000-1.0003.712厂2.0001.072实验目的掌握解线性方程组的Gauss消去法.基础理论3.000卜1]4.623x25.643人x3J[L000]=2.0000.000/Gauss消去法是通过对增广矩阵的初等行变换，将方程组变成上三角方程组，然后通过回代，从后到前依次求出各未知数.Gauss消去法的第k步（1MkMn-1）消元:若九丰0,则依次将增广矩阵第k行的-七/九倍加到第1行（k+1<i<n）,将第k列对角线下的元素都化成0.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程结果与分析附录：程序清单

《数值计算方法》实验3报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验3解线性方程组的直接法（之Gauss列主元消去法）实验题目用Gauss列主元消去法求解线性方程组'0.001-1.000厂2.0002.0003.7121.072'0.001-1.000厂2.0002.0003.7121.0723.000卜1]4.623x25.643人x3J[L000]=2.0000.000/掌握解线性方程组的Gauss列主元消去法.基础理论Gauss列主元消去法也是通过对增广矩阵的初等行变换，将方程组变成上三角方程组，然后通过回代，从后到前依次求出各未知数.Gauss列主元消去法的第k步（1MkMn-1）消元:先在九,a^,,ak中找绝对值最大的，将它所在的行与第k行交换，然后将第k行的-七/%倍加到第i行（k+1<i<n），将第k列对角线下的元素都化成0.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程结果与分析附录：程序清单

《数值计算方法》实验3报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验3解线性方程组的直接法（之Doolittle分解）实验题目对矩阵A进行Doolittle分解，其中(621-1、2410A=114-1[-10-13；3.实验目的掌握矩阵的Doolittle分解.4.基础理论矩阵的Doolittle分解是指将矩阵A=（七）呻可以分解为一个单位下三角矩阵和一个上三角矩阵的乘积.若设(100…0、'uuu…u、1112131nl10…00uu…u2122232nL=ll1…0,U=00u…u3132333n•••••.....・----.・...・・・.・12ln3…1^V000…u/nn则可依如下顺序公式计算1「九=akj-如l*r=1j=k,k+1,…，n1l=(aikik—Jlu)/u,ittkkkt=1i=k+1,k+2,…，n其中k=1,2,…，n.5.实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程（1）按计算公式依次计算一行u同时计算一列l；（2）因为计算完u..（或l..）后，a..就不再使用，为节省存储空间，将计算的u..（和l..）ijljljlj.j仍存放在矩阵A中的相应位置；（3）使用L矩阵和U矩阵时需要根据元素所在位置取固定值或A中相应位置的值.L对角线上的元素为1,上三角部分为0,下三角部分为A中对应的元素；U的下三角部分为0,上三角部分为A中对应的元素.

算法的伪代码描述：Input:阶数n；矩阵AOutput：矩阵L和U（合并存储在数组A中）行号伪代码注释1forkj1ton2forjjkton3u=a-*11u,kjkjkrrjr=1计算U的第k行（计算结果存放在aj4endfor5forijkton6l=(a-^11u)/uikikittkkkt=1计算L的第k列（计算结果存放在aj7endfor8endfor9returnL,U输出L和U结果与分析附录：程序清单rvf““““““““““““““““““““““““““““““““““““““““““““““““““““%**个*************************************************%程序名:Doolittle.m*%程序功能：矩阵LU分解中的Doolittle分解.rvf““““““““““““““““““““““““““““““““““““““““““““““““““““%****************************************************clc;%矩阵阶数clearall;n=4;A=[621-1;2410;114-1;-10-13]disp('A=');disp(A);%LU分解(Doolittle分解)fork=1:n%计算矩阵U的元素u_{kj}%（可参照下面l_{ik}的公式填写%矩阵阶数知i=k-1:nfcri=1:k-1Ag=A(ik)-A&n项均％计算矩阵L的元素l_(ik}endend%L在A下三角,U在上三角(对角线为1)enddisp('分解结果：');disp('L=');fori=1:nforj=1:nifi>j%在下三角部分，则取A对于的元素显示fprintf('%8.4f,A(i,j));elseifi==j%在对角线上，则显示1fprintf('%8d',1);else%在上三角部分，则显示0fprintf('%8d',0);endendfprintf('\n');%换行enddisp('U=');fori=1:nforj=1:nifi<=j%在上三角部分或对角线上，则取A对于的元素显示fprintf('%8.4f,A(i,j));else%在下三角部分，则显示0fprintf('%8d',0);endendfprintf('\n');%换行end《数值计算方法》实验3报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验3解线性方程组的直接法（之LU分解法）实验题目用LU分解（Doolittle分解）法求解线性方程组x+2x+6x=1<2x+5x+15x=36x'+15；+46【=10

、123实验目的熟悉解线性方程组LU分解法.基础理论、若将矩阵A进行了Doolittle分解，A=LU,则解方程组杰=b可以分解求解两个三角方程组Ly=b和Ux=M它们都可直接代入求解，其中Ly=b的代入公式为y=b一云ly,k=1,2,…,nj=1而ux=y的代入公式为x=（y-^Lux）/u,k=n,n-1,•…，1.kkkjjkkj=k+1实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程（1）Doolittle分解过程依次计算一行u同时计算一列l完成，并将计算的勺（和l〃）仍存放在矩阵A中的相应位置；""（2）求解方程组的代入公式中用到的勺和1〃都直接在A的相应位置取值即可.算法的伪代码描述：""Input：阶数n；矩阵A；常数项向量bOutput:方程组Ax=b的解x行号伪代码注释1LU分解参加上一实验程序2fork,1ton3y=b一支11y,k=1,2,…，nj=1代入法求解Ly=b4endfor

5fork,1ton6x=(y—才ux)/u,k=n,n-1,,1kkkjjkkj=k+1回代求解Ux=y7endfor8returnX],…,xn7.结果与分析8.附录：程序清单rvf$$$$$$$力e“$$$$$$$$$$$$$$““$$$$$$$$e%****************************************************%程序名：LinearSystemByLU.m*%程序功能：利用LU分解(Doolittle分解)解方程组.rvf$$$$$$-力、e“$$$$$$$$$$$$$$““$$$$$$$$e%****************************************************clc;%矩阵阶数clearall;n=3;A=[126;2515;61546];b=[1;3;10];%LU分解(Doolittle分解)fork=1:n%计算矩阵U的元素u_{kj}%（可参照下面l_{ik}的公式填写%矩阵阶数融i=k-1:nfcrj=1:k-1%计算矩阵L的元素l_{ik}AUk)=A(就-A(i:ij*A(j,k);endAMk)=AMk)'AfKk):%计算矩阵L的元素l_{ik}end%L在A下三角，U在上三角（对角线为1）endfork=1:ny(k)=b(k);%用代入法求解下三角方程组fork=1:ny(k)=b(k);%用代入法求解下三角方程组Ly=b%七=bk-11/.,k=1,2,…,nj=13enddisp('方程组Ly=b的解：y=');disp(y');fork=n:-1:1%回代求解上三角方程组Ux=y一.寸%Xk=y一乙x(k)=y(k);ux,k=n,n-1,•…,1j=k+1enddisp('原方程组的解：x=');一.寸%Xk=y一乙《数值计算方法》实验3报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验3解线性方程组的直接法(之Cholesky分解)实验题目对矩阵A进行Cholesky分解，其中f621-1]A=241°.114-1"-10-13J实验目的理解矩阵的Cholesky分解.基础理论矩阵的Cholesky分解是指将矩阵A=(%)心可以分解为一个下三角矩阵L和L转置的乘积，即A=LLt，其中L各元素可依如下顺序公式计算I七="_芸%＜》r=1l=(a-习ll)/1,i=k+1,k+2,,nikikittkkkIt=1其中k=1,2,...,n.实验环境操作系统：Windowsxp；程序设计语言：VC++实验过程按计算公式依次先计算一列对角线上的元素1kk，再计算这列其他元素履且对称位置的元素也取同一个值；因为计算完1..后，a„就不再使用，为节省存储空间，将计算的1..仍存放在矩阵A中的相应位置；““"使用L矩阵时需要根据元素所在位置取固定值或A中相应位置的值.L上三角部分为0，对角线和下三角部分为A中对应的元素.算法的伪代码描述：Input:阶数n；矩阵AOutput：矩阵L(合并存储在数组A中)行号伪代码注释1fork—1ton

2l=\a-t^112kk\|kkkr'r=13fori,kton41=(a-t111)/1ikikittkkkt=1计算结果存放在a..5endfor6endfor7returnL输出L结果与分析附录：程序清单%**个**********************************************************%*程序名：Cholesky.m*%*程序功能：对称正定矩阵的Cholesky分解.%*************************************************************n=4;%矩阵阶数A=[6,2,1,-1;2,4,1,0;1,1,4,-1;-1,0,-1,3];disp('A=');fori=1:nforj=1:nfprintf('%10.4f,A(i,j));%共占14位endfprintf('\n');%一行结束换行endfork=1:n%计算对角线上的Zfork=1:n%计算对角线上的Z_{kk}%计算其他的匚{ik}%和匚{ki}for户1:虹]A(k,k)=A-(l£nk)-A(Kjr2：endA'tki=，qrt[A[哄由血i=k+l:nfcrj=l:k-lA(i,k)=A(i,k)-AM+A03k)：endAftk)=Aftt)/A(ktXA(Ki)=A(i,k)：end

end%L在A下三角，Lend%L在A下三角，LAT在上三角forj=1:nifi>=j%在下三角部分或对角线上，则取A对于的元素显示fprintf('%10.4f,A(i,j));else%在上三角部分，则显示0fprintf('%10d',0);endendfprintf('\n');%换行end《数值计算方法》实验3报告班级：20xx级XXXXx班学号：20xx2409xxxx姓名：XXX成绩：实验名称实验3解线性方程组的直接法(之改进的Cholesky分解)实验题目对矩阵A进行改进的Cholesky分解，其中(621-1)2410A=114-1-1V10-13°7实验目的理解矩阵改进的Cholesky分解.基础理论矩阵的改进的Cholesky分解是指将矩阵A=(%)心可以分解为一个单位下三角矩阵L和对角矩阵D及L转置的乘积，即A=LDLt，其中L和D各元素可依如下顺序公式计算d=a一如d12kkkrkrr=11=(a-切dl1)/d,i=k+1,k+2,,nikiktitktkIt=1其中k=1，2，…，n.5.实验环境操作系统：Windowsxp；程序设计语言：VC++6.实验过程按计算公式依次先计算D的一个元素dk，再计算L中这列的元素匕/且对称位置的元素也取同一个值；因为计算完dk和/〃后，akk或％就不再使用，为节省存储空间，将计算的％或1..仍存放在矩阵A中的相应位置；使用L矩阵时需要根据元素所在位置取固定值或A中相应位置的值.L对角线和上三角部分为0，下三角部分为A中对应的元素；D对角线为A中对应的元素，其余都是0.算法的伪代码描述：Input:阶数n；矩阵AOutput：矩阵L(合并存储在数组A中)行号伪代码注释1forkj1ton

2d=a一切d12r=13fori,kton41=(a-云dl1)/dikiktitktkt=1计算结果存放在a--IJ5endfor6endfor7returnL,D输出L和D7.结果与分析8.附录：程序清单%**个**********************************************************%*程序名：ImprovedCholesky.m*%*程序功能：对称正定矩阵的改进的Cholesky分解.%*************************************************************n=4;%矩阵阶数A=[6,2,1,-1;2,4,1,0;1,1,4,-1;-1,0,-1,3];disp('A=');fori=1:nforj=1:nfprintf('%10.4f,A(i,j));%共占14位endfprintf('\n');%一行结束换行end%Cholesky分解fork=1:n%计算D对角线上的w_{kk}%计算L的元素Z_{ik}%和L转置的元素Z_{ki}%计算D对角线上的w_{kk}%计算L的元素Z_{ik}%和L转置的元素Z_{ki}endfori=k+l:nA(U)=A{U)-Ag*削ij)*endA(Lk)=A(iX)..0«k〕；A(ti)=A(U)；end

end%L在A下三角,Dend%L在A下三角,D在对角线forj=1:nifi>j%在下三角部分，则取A对于的元素显示fprintf('%10.4f,A(i,j));elseifi==j%在对角线上，则显示1fprintf('%10d',1);else%在上三角部分，则显示0fprintf('%10d',0);endendfprintf('\n');%换行enddisp('D=');fori=1:nforj=1:nifi==j%在对角线上，则取A对于的元素显示fprintf('%10.4f,A(i,j));else%其余显示0fprintf('%10d',0);endendfprintf('\n');%换行end《数值计算方法》实验3报告班级：20xx级XXXXx班姓名：XXX学号：20xx2409xxxx成绩：对于系数矩阵为三对角矩阵的方程组其Crout分解可分解为(b1a2C1b2a对于系数矩阵为三对角矩阵的方程组其Crout分解可分解为(b1a2C1b2an-1n-1ancn-1bnJ(s1a2asn-1n-1asnnt11t2实验3解线性方程组的直接法（之追赶法）实验题目用追赶法求解线性方程组-13-10实验目的熟悉解线性方程组的追赶法.基础理论这样，解方程组可以由如下2步完成:“追”：s=s=b,t=c/s,s=b—at,i-1i-1i-1iiii-1"/s1，J=(f—ay)/s,i=2,3,•••,n,iiii-1i其中：（匕，…，顷为方程组的常数项，八没用；“赶”：Xn=yn,x=y-1x,i=n-1,n-2,—,1.iiii+1实验环境操作系统：Windowsxp；程序设计语言：Matlab实验过程在“追”的过程中，向量s和y都有n个元素，t只有n-1个元素，又七和"的计算公式与其它s和y不同，所以先单独计算s和y，然后在一个n-1次循环中，求其它s和y以ii11ii

由于在“追”的过程中，b,c和f在分别计算完对应的s,t.和y后就不再使用，所iiiiIi以借用数组b，c和f存储向量s，t和y；同样在“赶”的过程中，七在计算完对应的X后就不再使用，所以再一次借用数组f存储向量x追赶法算法的伪代码描述：Input:阶数n；三对角矩阵的三条对角线向量a，b，c，常数项向量fOutput:方程组的解x行号伪代码注释1s1。b1；追2y】。f1/s1；3forkj2ton4tk-1jck-1/sk-1；5skjbk-akL；6ykj(fk-akyk-1)/sk；7endfor8