数值计算实验课题目
数值计算实验报告
本科实验报告课程名称：计算机数值方法实验项目：方程求根、线性方程组的直接解法、线性方程组的迭代解法、代数插值实验地点：专业班级：学生姓名：指导教师：实验一方程求根}五、实验结果与分析二分法实验结果迭代法实验结果结果分析：本题目求根区间为[1，2]，精度满足|x*-x n|<0.5×10-5，故二分法用公式|x*-x n|<(b-a)/ 2n，可求得二分次数并输出每次结果。
对迭代法首先要求建立迭代格式。
迭代格式经计算已输入程序之中，故直接给初值便可利用迭代法求出精度下的解。
六、讨论、心得每次的实验都是对已学过的理论知识的一种实战。
通过本次实验，我将二分法与迭代法的思路清晰化并且将其变成计算机设计语言编写出来，运用到了实际解决问题上感觉很好。
我自认为本次跟其他同学比较的优点在于我在二分法实现的时候首先利用换底公式将需要的二分次输输出，如此便很清晰明了的知道接下来每一步的意思。
迭代法给我的感觉便是高度的便捷简化，仅用几行代码便可以同样解决问题。
相比较二分法来说，我更喜欢迭代的思路。
实验二线性方程组的直接解法for(k=n-2;k>=0;k--){sum=0;for(j=k+1;j<n;j++)sum=sum+a[k][j]*x[j];x[k]=(b[k]-sum)/a[k][k];}for(i=0;i<n;i++)printf("x[%d]=%f ",i,x[i]); printf("\n"); //输出解向量x}五、实验结果与分析结果结果分析：如上图所示，输入线性方程组元数n=3，则会要求输入3*3的系数矩阵A与向量b构成的增广矩阵。
根据算法需要将系数矩阵A消元成上三角矩阵。
随后根据矩阵乘法公式变形做对应的回代。
六、讨论、心得本次实验在编写时候感觉还好，感觉将思路变成了程序设计语言，得以实现题目的要求。
但是在运行以及结果分析的时候，感觉到了本实验的一些不足之处：就是我的实验虽然可以实现不同的元数的线性方程组求解，但是缺少了分析初始条件——主元素不能为零。
数值计算方法课后习题答案
第一章 绪论（12）1、设0>x ，x 的相对误差为δ，求x ln 的误差。
[解]设0*>x 为x 的近似值，则有相对误差为δε=)(*x r ，绝对误差为**)(x x δε=，从而x ln 的误差为δδεε=='=*****1)()(ln )(ln x x x x x ， 相对误差为****ln ln )(ln )(ln x x x x rδεε==。
2、设x 的相对误差为2%，求n x 的相对误差。
[解]设*x 为x 的近似值，则有相对误差为%2)(*=x r ε，绝对误差为**%2)(x x =ε，从而nx 的误差为nn x x nxn x x n x x x **1***%2%2)()()()(ln *⋅=='=-=εε，相对误差为%2)()(ln )(ln ***n x x x nr==εε。
3、下列各数都是经过四舍五入得到的近似数，即误差不超过最后一位的半个单位，试指出它们是几位有效数字：1021.1*1=x ，031.0*2=x ，6.385*3=x ，430.56*4=x ，0.17*5⨯=x 。
[解]1021.1*1=x 有5位有效数字；0031.0*2=x 有2位有效数字；6.385*3=x 有4位有效数字；430.56*4=x 有5位有效数字；0.17*5⨯=x 有2位有效数字。
4、利用公式（3.3）求下列各近似值的误差限，其中*4*3*2*1,,,x x x x 均为第3题所给的数。
（1）*4*2*1x x x ++； [解]3334*4*2*11***4*2*1*1005.1102110211021)()()()()(----=⨯=⨯+⨯+⨯=++=⎪⎪⎭⎫ ⎝⎛∂∂=++∑x x x x x f x x x e nk k k εεεε；（2）*3*2*1x x x ；[解]52130996425.010********.2131001708255.01048488.2121059768.01021)031.01021.1(1021)6.3851021.1(1021)6.385031.0()()()()()()()()(3333334*3*2*1*2*3*1*1*3*21***3*2*1*=⨯=⨯+⨯+⨯=⨯⨯+⨯⨯+⨯⨯=++=⎪⎪⎭⎫⎝⎛∂∂=-------=∑x x x x x x x x x x x f x x x e n k k kεεεε；（3）*4*2/x x 。
丁丽娟《数值计算方法》五章课后实验题答案（源程序很详细,且运行无误）
丁丽娟《数值计算方法》五章课后实验题答案（源程序都是自己写的，很详细，且保证运行无误）我做的五章数值实验作业题目如下：第二章：1、2、3、4题第三章：1、2题第四章：1、2题第六章：2、3题第八章：1、2题第二章1：(1) 对A进行列主元素三角分解：function [l u]=myfun(A) n=size(A); for k=1:n for i=k:n sum=0; m=k; for j=1:(k-1) sum=sum+A(i,j)*A(j,k); end s(i)=A(i,k)-sum; if abs(s(m))<abs(s(i)) m=i; end end for j=1:n c=A(m,j); A(m,j)=A(k,j); A(k,j)=c; end for j=k:n sum=0; for r=1:(k-1) sum=sum+A(k,r)*A(r,j); end u(k,j)=A(k,j)-sum; A(k,j)=u(k,j); end for i=1:n l(i,i)=1; end for i=(k+1):n sum=0; for r=1:(k-1) sum=sum+A(i,r)*u(r,k); end l(i,k)=(A(i,k)-sum)/u(k,k); A(i,k)=l(i,k); end end 的列主元素三角分解：求A的列主元素三角分解：>>A=[1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; >>[L,U]=myfun(A) 结果：L = 1.0000 0 0 0 0 1.0000 1.0000 0 0 0 1.0000 0.5000 1.0000 0 0 1.0000 0.7500 0.7500 1.0000 0 1.0000 0.2500 0.7500 -1.0000 1.0000 U = 1.0000 1.0000 1.0000 1.0000 1.0000 0 4.0000 14.0000 34.0000 69.0000 0 0 -2.0000 -8.0000 -20.5000 0 0 0 -0.5000 -2.3750 0 0 0 0 -0.2500 (2) 求矩阵的逆矩阵A -1: inv(A) 结果为：ans = 5 -10 10 -5 1 -10 30 -35 19 -4 10 -35 46 -27 6 -5 19 -27 17 -4 1 -4 6 -4 1 （3）检验结果：E=diag([1 1 1 1 1]) A\E ans = 5 -10 10 -5 1 -10 30 -35 19 -4 10 -35 46 -27 6 -5 19 -27 17 -4 1 -4 6 -4 1 2: 程序：程序：function d=myfun(a,b,c,d,n) for i=2:n l(i)=a(i)/b(i-1); a(i)=l(i); u(i)=b(i)-c(i-1)*a(i); b(i)=u(i); y(i)=d(i)-a(i)*d(i-1); d(i)=y(i); end x(n)=d(n)/b(n); d(n)=x(n); for i=(n-1):-1:1 x(i)=(d(i)-c(i)*d(i+1))/b(i); d(i)=x(i); end 求各段电流量程序：求各段电流量程序：for i=2:8 a(i)=-2; end b=[2 5 5 5 5 5 5 5]; c=[-2 -2 -2 -2 -2 -2 -2]; V=220; R=27; d=[V/R 0 0 0 0 0 0 0]; n=8; I=myfun(a,b,c,d,n) 运行程序得：运行程序得：I = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477 3：程序：（1）求矩阵A和向量b的matlab程序：function [A b]=myfun(n) for i=1:n X(i)=1+0.1*i; end for i=1:n for j=1:n A(i,j)=X(i)^(j-1); end end for i=1:n b(i)=sum(A(i,:)); end 求n=5时A1，b1及A1的2-条件数程序运行结果如下：条件数程序运行结果如下： n=5；[A1,b1]=myfun(n) A1 = 1.0000 1.1000 1.2100 1.3310 1.4641 1.0000 1.2000 1.4400 1.7280 2.0736 1.0000 1.3000 1.6900 2.1970 2.8561 1.0000 1.4000 1.9600 2.7440 3.8416 1.0000 1.5000 2.2500 3.3750 5.0625 b1 = 6.1051 7.4416 9.0431 10.9456 13.1875 cond2=cond(A1,2)cond2 = 5.3615e+005 条件数程序运行结果如下：求n=10时A2，b2及A2的2-条件数程序运行结果如下：n=10; [A2,b2]=myfun(n) A2 = 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.7716 1.9487 2.1436 2.3579 1.0000 1.2000 1.4400 1.7280 2.0736 2.4883 2.9860 3.5832 4.2998 5.1598 1.0000 1.3000 1.6900 2.1970 2.8561 3.7129 4.8268 6.2749 8.1573 10.6045 1.0000 1.4000 1.9600 2.7440 3.8416 5.3782 7.5295 10.5414 14.7579 20.6610 1.0000 1.5000 2.2500 3.3750 5.0625 7.5938 11.3906 17.0859 25.6289 38.4434 1.0000 1.6000 2.5600 4.0960 6.5536 10.4858 16.7772 26.8435 42.9497 68.7195 1.0000 1.7000 2.8900 4.9130 8.3521 14.1986 24.1376 41.0339 69.7576 118.5879 1.0000 1.8000 3.2400 5.8320 10.4976 18.8957 34.0122 61.2220 110.1996 198.3593 1.0000 1.9000 3.6100 6.8590 13.0321 24.7610 47.0459 89.3872 169.8356 322.6877 1.0000 2.0000 4.0000 8.0000 16.0000 32.0000 64.0000 128.0000 256.0000 512.0000 b2 = 1.0e+003 * 0.0159 0.0260 0.0426 0.0698 0.1133 0.1816 0.2866 0.4451 0.6801 1.0230 cond2=cond(A2,2) cond2 = 8.6823e+011 条件数程序运行结果如下：求n=20时A3，b3及A3的2-条件数程序运行结果如下：n=20; [A3,b3]=myfun(n) A3 = 1.0e+009 * Columns 1 through 10 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 11 through 20 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0003 0.0005 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0003 0.0006 0.0013 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0003 0.0007 0.0015 0.0032 0.0000 0.0000 0.0000 0.0001 0.0001 0.0003 0.0006 0.0014 0.0032 0.0075 0.0000 0.0000 0.0000 0.0001 0.0002 0.0005 0.0012 0.0029 0.0070 0.0167 0.0000 0.0000 0.0001 0.0001 0.0004 0.0009 0.0023 0.0058 0.0146 0.0364 0.0000 0.0000 0.0001 0.0002 0.0006 0.0017 0.0044 0.0113 0.0295 0.0766 0.0000 0.0001 0.0002 0.0004 0.0011 0.0030 0.0080 0.0215 0.0581 0.1570 0.0000 0.0001 0.0002 0.0007 0.0018 0.0051 0.0143 0.0400 0.1119 0.3133 0.0000 0.0001 0.0004 0.0010 0.0030 0.0086 0.0250 0.0726 0.2105 0.6103 0.0001 0.0002 0.0005 0.0016 0.0048 0.0143 0.0430 0.1291 0.3874 1.1623 b3 = 1.0e+009 * Columns 1 through 10 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0002 0.0004 0.0010Columns 11 through 20 0.0025 0.0059 0.0132 0.0287 0.0606 0.1246 0.2494 0.4874 0.9316 1.7434 cond2=cond(A3,2) cond2 =3.2395e+022 由上述运行结果可知：它们是病态的，而且随着n的增大，矩阵的病态变得严重。
实验五+MATLAB数值计算(含实验报告)
实验五 MATLAB 数值计算一、实验目的1．掌握求数值导数和数值积分的方法。
2．掌握代数方程数值求解的方法。
3．掌握常微分方程数值求解的方法。
二、实验的设备及条件计算机一台（带有MATLAB7.0以上的软件环境）。
设计提示1．参考本节主要内容，学习并理解相关函数的含义及调用方法。
三、实验内容1．线性系统方程：分别使用左除（\）和求逆（inv ）求解下面系统方程的解：⎪⎩⎪⎨⎧=+=+=++377251463c b b a c b a2. 数值积分：使用quad 和trapz 求解⎰-503/dx xe x 的数值积分，并与其解析解9243/5+--e 相比较；3. 请完成教材P154页中实验指导环节的实验内容第2题4. 请完成教材P155页中思考练习的第3题（1），并绘制解在该求解区间（即[0,5]）上的图像；。
5、请完成教材P164页实验指导环节的实验内容第5题。
（提示：该函数的符号导数，可以通过函数diff 求得。
首先定义符号变表达式，如求sin(x)的一阶符号导数，可以先定义f=’sin(x)’;df=diff(f);可求得df=cos(x)。
其中df 即为函数f 的一阶符号导数）。
四、实验报告要求（包含预习报告要求和最终报告要求）1.实验名称2.实验目的3.实验设备及条件4.实验内容及要求5.实验程序设计指程序代码。
6.实验结果及结果分析实验结果要求必须客观,现象。
结果分析是对实验结果的理论评判。
7.实验中出现的问题及解决方法8. 思考题的回答五、实验报告的提交方式Word文档，命名方式:实验号_你的学号_姓名例如本次实验：实验一_000000001_张三.doc（信息101提交报告邮箱）：E_mail: *******************(网络工程101提交作业邮箱)：E_mail: *******************（注意网络班的M是大写的)下一次课前提交,过期不收!六、参考文献参考教材和Matlab帮助文件。
数值计算课后习题答案
习 题 一 解 答1．取3.14，3.15，227，355113作为π的近似值，求各自的绝对误差，相对误差和有效数字的位数。
分析：求绝对误差的方法是按定义直接计算。
求相对误差的一般方法是先求出绝对误差再按定义式计算。
注意，不应先求相对误差再求绝对误差。
有效数字位数可以根据定义来求，即先由绝对误差确定近似数的绝对误差不超过那一位的半个单位，再确定有效数的末位是哪一位，进一步确定有效数字和有效数位。
有了定理2后，可以根据定理2更规范地解答。
根据定理2，首先要将数值转化为科学记数形式，然后解答。
解：（1）绝对误差:e(x)=π－3.14＝3.14159265…－3.14＝0.00159…≈0.0016。
相对误差：3()0.0016()0.51103.14r e x e x x-==≈⨯有效数字：因为π＝3.14159265…=0.314159265…×10，3.14＝0.314×10，m=1。
而π－3.14＝3.14159265…－3.14＝0.00159…所以│π－3.14│＝0.00159…≤0.005=0.5×10－2＝21311101022--⨯=⨯所以，3.14作为π的近似值有3个有效数字。
（2）绝对误差:e(x)=π－3.15＝3.14159265…－3.14＝－0.008407…≈－0.0085。
相对误差：2()0.0085()0.27103.15r e x e x x--==≈-⨯有效数字：因为π＝3.14159265…=0.314159265…×10，3.15＝0.315×10，m=1。
而π－3.15＝3.14159265…－3.15＝－0.008407…所以│π－3.15│＝0.008407……≤0.05=0.5×10－1＝11211101022--⨯=⨯所以，3.15作为π的近似值有2个有效数字。
（3）绝对误差:22() 3.14159265 3.1428571430.0012644930.00137e x π=-=-=-≈-相对误差：3()0.0013()0.4110227r e x e x x--==≈-⨯有效数字：因为π＝3.14159265...=0.314159265 (10)22 3.1428571430.3142857143107==⨯，m=1。
数值计算方法课后实验-李维国
实验3.1 Gauss消去法的数值稳定性实验实验目的：观察和理解高斯消元过程中出现小主元(即|a(k)kk|很小)时引起方程组解数值不稳定性．实验内容：求解方程组Ax=b,其中实验要求：（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的．（2）用高斯列主元消去法求得L和U及解向量x1,x2 ．（3）用不选主元的高斯消去法求得L和U及解向量x1,x2 ．（4）观察小主元并分析对计算结果的影响．解：（1）cond(A1)=68.4296cond(A2)=8.9939A1矩阵条件数远大于1，A1矩阵为病态;A2矩阵条件数没有远远大于1，A2矩阵为良态。
（2）高斯列主元消去法：程序如下ClearA1=[0.3*power(10,-15),59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];b1=[59.17;46.18;1;2];% A1=[1,2,-1,1;1,1,2,-1;3,-1,1,1;2,1,3,-1];%b1=[1;-2;6;-1];A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2];b2=[8;5.900000000001;5;1];A=A1; b=b1;n=input('n=');for k=1:n-1 [a,b3]=max(A(k:n,k)); A([k k-1+b3],:)=A([k-1+b3 k],:); b([k k-1+b3],:)=b([k-1+b3 k],:); if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)./A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); else stop end endA; b; L=tril(A,-1);U=triu(A,0); for i=1:1:n L(i,i)=1;end L; for j=1:n-1 b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);end b(n)=b(n)/L(n,n); b; y=b;for j=n:-1:2 y(j)=y(j)/U(j,j); y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);end y(1)=y(1)/U(1,1);结果如下:方程组一：L1： U1：解向量x1结果如下：方程组二结果如下：L2：U2:解向量x2结果如下：（3）不选主元的分解程序如下：ClearA1=[0.3*power(10,-15),59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1]; b1=[59.17;46.18;1;2];%A1=[6,2,1,-1;2,4,1,0;1,1,4,-1;-1,0,-1,3];%b1=[6;-1;5;-5];A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2];b2=[ 8;5.900000000001;5;1]; A=A2; b=b2; n=input('n='); for k=1:n-1A(k+1:n,k)= A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); end L=tril(A,-1); U=triu(A,0); for i=1:1:n L(i,i)=1; end L ;for j=1:n-1 b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j); end b(n)=b(n)/L(n,n); b; y=b;for j=n:-1:2 y(j)=y(j)/U(j,j); y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j); end y(1)=y(1)/U(1,1);方程组一结果：L1: U1:解向量x1结果如下：方程组二结果：L2: U2:解向量x2结果如下：(4) 观察方程在两种不同方法下的结果可知：由于计算机字长是一定的，小主元会造成大数除以小数的结果超出字长，结果发生很大的变化。
数值计算大作业
课程设计课程名称：设计题目：学号：姓名：完成时间：题目一：非线性方程求根 一 摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。
本实验通过使用常用的求解方法二分法和Newton 法及改进的Newton 法处理几个题目，分析并总结不同方法处理问题的优缺点。
观察迭代次数，收敛速度及初值选取对迭代的影响。
用Newton 法计算下列方程（1） 310x x --= ， 初值分别为01x =，00.45x =，00.65x =； （2） 32943892940x x x +-+= 其三个根分别为1,3,98-。
当选择初值02x =时给出结果并分析现象,当6510ε-=⨯，迭代停止。
解：1)采用MATLAB 进行计算；首先定义了Newton 法：function kk=newton(f,df,x0,tol,N)% Newton Method （牛顿法）% The first parameter f is a external function with respect to viable x.（第一个参数也就是本题所用的函数f ）% The second parameter df is the first order diffential function of fx.（第二个参数也就是本体所用函数f 的导数方程df ） % x0 is initial iteration point(初值). % tol is the tolerance of the loop （精度）.% N is the maximum number of iterations （循环上限）. x=x0;f0=eval(f);df0=eval(df); n=0;disp(' [ n xn xn+1 fn+1 ]'); while n<=N x1=x0-f0/df0; x=x1; f1=eval(f); X=[n,x0,x1,f1]; disp(X);if abs(x0-x1)<tolfprintf('The procedure was successful.') kk=X; return else n=n+1; x0=x1;f0=f1;endendif n==N+1fprintf('the method failed after N iterations. '),kk=0;End我们把Newton法存为.m格式的文件；之后我们运行程序：clear;clc;syms xf=x^3-x-1;df=diff(f,x);x=newton(f,df,1,0.0001,50);x会得到一下结果[ n xn xn+1 fn+1 ]0 1.0000 1.5000 0.87501.0000 1.5000 1.0625 -0.86302.0000 1.0625 1.4940 0.8408到第50次迭代时候会出现该问题：47.0000 1.4898 1.0814 -0.816748.0000 1.0814 1.4898 0.816749.0000 1.4898 1.0814 -0.816750.0000 1.0814 1.4898 0.8167the method failed after N iterations.x =0；同样测试x0=0.45、0.65得不出结果，判断出初值离真值太远，所以我们采用牛顿下山法进行计算迭代：我们定义了其中的f函数和df函数，并且分别存为.m格式的文件，其代码如下：f:function y=f(x)y=x^3-x-1;df:function y=df(x)y=3*x^2-1;之后我们定义newton下山法同时也存为.m的程序：function [x,i]=downnewton(f,df,x0,tol)k=0;i=1;disp(' [ n xn xn+1 fn+1 ]'); while(k==0)fx=feval('f',x0);dfx=feval('df',x0);t=0;u=1;while(t==0)dx=-fx/dfx;x1=x0+u*dx;fx1=feval('f',x1);fx0=feval('f',x0);if(abs(fx1)>abs(fx0));u=u/2;elset=1;endendX=[i,x0,x1,fx1];disp(X);if(abs(fx1)<tol)k=1;elsex0=x1;i=i+1;endendx=x1;i=i;end之后带入x0=0.45;downnewton('f','df',0.45,10^(-6))[ n xn xn+1 fn+1 ]1.0000 0.4500 -0.4155 -0.65622.0000 -0.4155 -0.5857 -0.61523.0000 -0.5857 -0.5754 -0.61514.0000 -0.5754 -0.5782 -0.61515.0000 -0.5782 -0.5773 -0.61516.0000 -0.5773 -0.5774 -0.61517.0000 -0.5774 -0.5773 -0.61518.0000 -0.5773 -0.5774 -0.61519.0000 -0.5774 -0.5774 -0.615110.0000 -0.5774 -0.5774 -0.615111.0000 -0.5774 1.3131 -0.049012.0000 1.3131 1.3248 0.000513.0000 1.3248 1.3247 0.0000ans =1.3247带入x0=0.6;downnewton('f','df',0.6,10^(-6))[ n xn xn+1 fn+1 ]1.0000 0.6000 1.1406 -0.65662.0000 1.1406 1.3668 0.18663.0000 1.3668 1.3263 0.00674.0000 1.3263 1.3247 0.00005.0000 1.3247 1.3247 0.0000ans =1.3247带入x0=1;downnewton('f','df',1,10^(-6))[ n xn xn+1 fn+1 ]1.0000 1.0000 1.5000 0.87502.0000 1.5000 1.3478 0.10073.0000 1.3478 1.3252 0.00214.0000 1.3252 1.3247 0.0000ans =1.32472)同样采用Newton下山法：重新定义f、df:f:function y=f(x)y=x^3+94*x^2-389*x+294;df:function y=df(x)y=3*x^2+188*x-389;再带入初值x0=2;downnewton('f','df',2,5*10^(-6))[ n xn xn+1 fn+1 ]1 2 -98 0ans =-98得出x=-98；分析：先画出该函数的图像；x=(-100:.1:100);ezplot('x^3+94*x^2-389*x+294',[-100 100]) 得出该图像如图：-100-80-60-40-20020406080024681012141618x 105xx 3+94 x 2-389 x+294根据牛顿法的几何解释，在x0=2的点做切线，与y 相交，交点的横坐标值为x=-98则结束了该现象。
数值计算方法实验1
学院（系）名称：)()()()(0101112x x x f x f x f x x ---=附录（源程序及运行结果）：一．二分法#include<stdio.h>#include<math.h>double f(double x){return x*x-x-1;}void main(){float a=0,b=0,x=1,m,e;int k;while(f(a)*f(b)>0){printf("请输入区间a,b的值。
以及精度e\n");scanf("%f,%f,%f",&a,&b,&e);}k=0;if(f(a)*f(b)==0){if(f(a)==0)printf("使用二分法输出：a=%f,k=%d\n",a,k);elseprintf("使用二分法输出：b=%f,k=%d\n",b,k);}else{while(f(a)*f(b)!=0){m=(a+b)/2;if(fabs(a-b)/2<e){printf("使用二分法输出：m=%f,k=%d\n",m,k);break;}else {if(f(a)*f(m)>0)a=m;else b=m;k=k+1;}}}}运行结果：二．迭代法与牛顿迭代法#include<stdio.h>#include<math.h>double f(double x){return exp(-x);}double f1(double x){return (x*exp(x)-1);}double ff(double x){return (exp(x)+x*exp(x));}void diedaifa(double x0,double e,int N){double x1;int k=1;while(k!=N){x1=f(x0);if(fabs(x1-x0)>=e){k++;if(k==N)printf("迭代失败！\n");x0=x1;}else{printf("使用迭代法输出结果：%lf\n",x1);break;}}}void NDdiedaifa(double x0,double e,int N){int k=1;double x1;while(k!=N){if(ff(x0)==0)printf("公式f(x)奇异！\n");else{x1=x0-f1(x0)/ff(x0);if(fabs(x1-x0)>=e){k++;if(k==N)printf("迭代失败！\n");x0=x1;}else{printf("使用牛顿迭代法输出结果：%lf\n",x1);break;}}}}void main(){double x0,e;int N;printf("请输入初值：");scanf("%lf",&x0);printf("精度：");scanf("%lf",&e);printf("以及判定迭代失败的最大次数N：");scanf("%d",&N);diedaifa(x0,e,N);NDdiedaifa(x0,e,N);}运行结果：四．双点弦截法#include<stdio.h>#include<math.h>double f(double x){return (x*x*x+3*x*x-x-9);}void main(){double x0,x1,x2,e;int N;int k=1;printf("请输入初值x0和x1：");scanf("%lf,%lf",&x0,&x1);printf("精度：");scanf("%lf",&e);printf("以及判定迭代失败的最大次数N：");scanf("%d",&N);while(k!=N){x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0));if(fabs(f(x2))>=e){k++;if(k==N)printf("迭代失败！\n");x0=x1;x1=x2;}else{printf("使用双点弦截法输出结果：%lf\n",x2);break;}}}运行结果：。
数值计算方法实验报告(含所有)
本科实验报告课程名称：计算机数值方法实验项目：计算机数值方法实验实验地点：虎峪校区致远楼B401专业班级：软件学院1217班学号：******xxxx 学生姓名：xxx指导教师：xxx2014 年 5 月21 日太原理工大学学生实验报告五、实验结果与分析二分法割线法分析：由程序知，使用二分法和割线法均能计算出方程的根，但利用割线法要比二分法计算的次数少，并且能够较早的达到精度要求。
相比之下，割线法程序代码量较少，精简明了。
六、讨论、心得本次数值计算方法程序设计实验从习题练习中跳脱出来，直接面对实用性较强的程序代码编写。
效果很好，不仅加深对二分法、割线法的理解，还加强了实际用运能力。
将理论知识成功地转化成实践结果。
实验地点虎峪校区致远楼B401指导教师xx太原理工大学学生实验报告l[i][k]=a[i][k];for(r=1;r<k;++r){l[i][k]-=l[i][r]*u[r][k];}l[i][k]/= u[k][k];}l[k][k]=1.0;}for(i=1;i<=n;++i){y[i] = b[i];for(j=1;j<i;++j){y[i]-=l[i][j]*y[j];}}for(i=n;i>0;--i){x[i] = y[i];for(j=i+1;j<=n;++j){x[i]-=u[i][j]*x[j];}x[i]/= u[i][i];}for(i=1;i<=n;++i){printf("%0.2lf\n",x[i]);}return 0;}五、实验结果与分析完全主元素消元法：列主元素消元法：LU分解法：分析：对于两种高斯解方程，完全主元素跟列主元素都是先消元、再回代，由程序段可以发现，始终消去对角线下方的元素。
即，为了节约内存及时效，可以不必计算出主元素下方数据。
列主元素消元法的算法设计上优于完全主元素消元法，它只需依次按列选主元素然后换行使之变到主元素位置，再进行消元即可。
实验09数值微积分与方程数值解(第6章)
实验09数值微积分与方程数值解（第6章）《数学软件》课内实验王平（第6章MATLAB数值计算）一、实验目的1.掌握求数值导数和数值积分的方法。
2.掌握代数方程数值求解的方法。
3.掌握常微分方程数值求解的方法。
二、实验内容1.求函数在指定点的数值导数某f(某)1程序及运行结果：某2某36某2某3某2,某1,2,3 022.用数值方法求定积分(1)I120cot24in(2t)21dt的近似值。
程序及运行结果：(2)I220ln(1某)d某1某2程序及运行结果：3.分别用3种不同的数值方法解线性方程组6某5y2z5u49某y4zu133某4y2z2u13某9y2u11程序及运行结果：4.求非齐次线性方程组的通解2某17某23某3某463某15某22某32某449某4某某7某22341程序及运行结果（提示：要用教材中的函数程序line_olution）：5.求代数方程的数值解(1)3某+in某-e某=0在某0=1.5附近的根。
程序及运行结果（提示：要用教材中的函数程序line_olution）：(2)在给定的初值某0=1，y0=1，z0=1下，求方程组的数值解。
in某y2lnz70y33某2z10某yz50程序及运行结果：26.求函数在指定区间的极值某3co某某log某(1)f(某)在(0,1)内的最小值。
e某程序及运行结果：332(2)f(某1,某2)2某1在[0,0]附近的最小值点和最小值。
4某1某210某1某2某2程序及运行结果：7.求微分方程的数值解，并绘制解的曲线某d2ydy5y0d某2d某y(0)0y'(0)0程序及运行结果（注意：参数中不能取0，用足够小的正数代替）：令y2=y,y1=y'，将二阶方程转化为一阶方程组：1'5yy1某1某y2'y2y1y(0)0,y(0)0218.求微分方程组的数值解，并绘制解的曲线y'1y2y3y'yy213y'0.51yy123y1(0)0,y2(0)1,y3(0)1程序及运行结果：3三、实验提示四、教程：第6章MATLAB数值计算（2/2）6.2数值微积分p1556.2.1数值微分1.数值差分与差商对任意函数f(某)，假设h>0。
数值计算(二分法、简单迭代法、Newton迭代法、弦截法(割线法、双点弦法))
-1.50508
12
-1.50504
3
-1.50625
8
-1.50449
13
-1.50506
4
-1.49688
9
-1.50479
14
-1.50507
表1-1
区间[-1.2,-0.9]
k
xk
k
xk
k
xk
0
-1.05
5
-0.998437
10
-1.00005
1
-0.975
6
-1.00078
11
-0.999976
13
1.69015
20
1.69028
7
1.68753
14
1.6902
表2-3
牛顿
初值-1.5结果x=-1.50507
k
xk
k
xk
1
-1.5
4
-1.50504
2
-1.50471
5
-1.50506
3
-1.50497
6
-1.50507
表3-1
初值-1结果x=-1.50507
k
x
1
-1
2
-1
表3-2
初值1.6结果x=1.69028
步骤：1.计算原函数的导数f’(x);构造牛顿迭代公式
2.计算,若f’(x0)=0,退出计算，否则继续向下迭代。
3.若|x1-x0|满足精度要求，x1即为方程的近似解。
2.4弦截法
思想:为加速收敛，改用两个端点都在变动的弦，用差商替代牛顿迭代公式的导数f’(x)。
步骤：1.构造双点弦法的公式
2.计算x2=x1-f(x1)(x1-x0)/f(x1)-f(x0);
数值计算方法丁丽娟课后习题答案
数值计算方法丁丽娟课后习题答案【篇一：北京理工大学数值计算方法大作业数值实验1】）书p14/4分别将区间[?10,10]分为100,200,400等份，利用mesh或surf命令画出二元函数的三维图形。
z=|??|+ ??+?? +??++??【matlab求解】[x,y]=meshgrid(-10:0.1:10);a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);[x,y]=meshgrid(-10:0.05:10);a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);[x,y]=meshgrid(-10:0.025:10); a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);（二）书p7/1.3.2数值计算的稳定性（i）取= ??c语言程序—不稳定解 +=ln1.2，按公式=?? （n=1,2，…） #includestdio.h#includeconio.h#includemath.hvoid main(){float m=log(6.0)-log(5.0),n;int i;i=1;printf(y[0]=%-20f,m); while(i20){n=1/i-5*m;printf(y[%d]=%-20f,i,n);m=n;i++;if (i%3==0) printf(\n); }getch();}（ii） c语言程序—稳定解≈??[ ??+?? +?? ??+??按公式 =??(??)#includestdio.h#includeconio.h#includemath.hvoid main(){float m=(1/105.0+1/126.0)/2,n; k=n,n-1,n-2，…）（【篇二：北京理工大学数值计算方法大作业数值实验4】 p260/1考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为= ?????????????????? ?? ??????????= ?????????????? ??曲线关于原点对称，取a=1，参数s的变化范围[-5,5]，容许误差限分别是，,和。
《数值分析》课程设计—16题
《数值分析》课程设计—作业实验一1.1水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。
由于旅途的颠簸，大家都很疲惫，很快就入睡了。
第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。
解：一、问题分析：对于本题，比较简单，我们只需要判断原来椰子的个数及每个人私藏了一份之后剩下的是否能被5除余1，直到最后分完。
对于第一个程序，n取2000；对于第二个程序，n取20001，就能得到我们想要的结果，即原先一共有15621个椰子，最终平均每人得4092个椰子。
1.2 当0,1,2,,100n =时，选择稳定的算法计算积分10d 10nxn xe I x e --=+⎰. 解：一、问题分析：由10d 10nxn xe I x e --=+⎰知： 1101001==+⎰dx I I 以及： )1(11010101010)1(1nnx x nx x n n n e ndx e dx e e e I I ----+-+-==++=+⎰⎰ 得递推关系： ⎪⎩⎪⎨⎧--=-=-+n nn I e n I I I 10)1(1101101， 但是通过仔细观察就能知道上述递推公式每一步都将误差放大十倍，即使初始误差很小，但是误差的传播会逐步扩大，也就是说用它构造的算法是不稳定的，因此我们改进上述递推公式（算法）如下：⎪⎪⎩⎪⎪⎨⎧--=-=+-))1(1(101)1(101110n n n I e n I I I通过比较不难得出该误差是逐步缩小的，即算法是稳定的。
二、问题求解：为了利用上面稳定的算法，需要我们估计初值100I 的值。
数值分析上机作业1-1解析
数值计算方法上机题目11、实验1. 病态问题实验目的：算法有“优”与“劣”之分，问题也有“好”和“坏”之别。
所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。
希望读者通过本实验对此有一个初步的体会。
数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。
病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。
问题提出：考虑一个高次的代数多项式∏=-=---=201)()20)...(2)(1()(k k x x x x x p （E1-1）显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。
现考虑该多项式方程的一个扰动0)(19=+xx p ε (E1-2)其中ε是一个非常小的数。
这相当于是对（E1-1）中19x 的系数作一个小的扰动。
我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。
实验内容：为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数u ＝roots （a ）其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。
设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程0...1121=++++-n n n n a x a x a x a的全部根，而函数b=poly(v)的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。
可见“roots ”和“Poly ”是两个互逆的运算函数.ve=zeros(1,21); ve(2)=ess;roots(poly(1:20))+ve)上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。
实验要求：（1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。
数值计算方法实验指导(Matlab版)
《数值计算方法》实验指导（Matlab版）学院数学与统计学学院计算方法课程组《数值计算方法》实验1报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验1 算法设计原则验证（之相近数相减、大数吃小数和简化计算步骤） 2. 实验题目(1) 取1610=z ，计算z z -+1和)1/(1z z ++，验证两个相近的数相减会造成有效数字的损失．(2) 按不同顺序求一个较大的数（123）与1000个较小的数（15310-⨯）的和，验证大数吃小数的现象．(3) 分别用直接法和九韶算法计算多项式n n n n a x a x a x a x P ++++=--1110)(在x =1.00037处的值．验证简化计算步骤能减少运算时间．对于第(3)题中的多项式P (x )，直接逐项计算需要2112)1(+=+++-+n n n 次乘法和n 次加法，使用九韶算法n n a x a x a x a x a x P ++++=-)))((()(1210则只需要n 次乘法和n 次加法． 3. 实验目的验证数值算法需遵循的若干规则． 4. 基础理论设计数值算法时，应避免两个相近的数相减、防止大数吃小数、简化计算步骤减少运算次数以减少运算时间并降低舍入误差的积累．两相近的数相减会损失有效数字的个数，用一个大数依次加小数，小数会被大数吃掉，乘法运算次数太多会增加运算时间． 5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab6. 实验过程(1) 直接计算并比较；(2) 法1：大数逐个加1000个小数，法2：先把1000个小数相加再与大数加； (3) 将由高次项到低次项的系数保存到数组A[n]中，其中n 为多项式次数．7. 结果与分析 (1) 计算的z z -+1= ，)1/(1z z ++．分析：(2) 123逐次加1000个6310-⨯的和是 ，先将1000个6310-⨯相加，再用这个和与123相加得．分析：(3) 计算次的多项式：直接计算的结果是，用时；用九韶算法计算的结果是，用时．分析：8. 附录：程序清单(1) 两个相近的数相减．%*************************************************************%* 程序名：ex1_1.m *%* 程序功能：验证两个相近的数相减会损失有效数字个数 *%*************************************************************z=1e16;x,y======================================================================(2) 大数吃小数%*************************************************************%* 程序名：ex1_2.m *%* 程序功能：验证大数吃小数的现象. *%*************************************************************clc; % 清屏clear all; % 释放所有存变量format long; % 按双精度显示浮点数z=123; % 大数t=3e-15; % 小数x=z; % 大数依次加小数% 重复1000次给x中加上ty=0; % 先累加小数% 重复1000次给y中加上ty=z + y; % 再加到大数x,y======================================================================(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,6 5,55,80,78,77,98,56];A(10001)=0; % 扩展到10001项，后面的都是分量0% A为多项式系数,从高次项到低次项x=1.00037;n=9000; % n为多项式次数% 直接计算begintime=clock; % 开始执行的时间 % 求x的i次幂% 累加多项式的i次项endtime=clock; % 完毕执行的时间time1=etime(endtime,begintime); % 运行时间disp('直接计算');disp(['p(',num2str(x),')=',num2str(p)]);disp([' 运行时间: ',num2str(time1),'秒']);% 九韶算法计算begintime=clock; % 开始执行的时间% 累加九韶算法中的一项endtime=clock; % 完毕执行的时间time2=etime(endtime,begintime); % 运行时间disp(' ');disp('九韶算法计算');disp(['p(',num2str(x),')=',num2str(p)]);disp([' 运行时间: ',num2str(time2),'秒']);《数值计算方法》实验1报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验1 算法设计原则验证（之数值稳定性） 2. 实验题目 计算定积分⎰==-1110,1,0,d n x e xI x nn ，分别用教材例1-7推导出的算法A 和B ，其中：算法A ：⎩⎨⎧≈-=-6321.0101I nI I n n 算法B ：⎪⎩⎪⎨⎧≈-=-0)1(1101I I nI n n 验证算法不稳定时误差会扩大．3. 实验目的验证数值算法需遵循的若干规则． 4. 基础理论设计数值算法时，应采用数值稳定性好的算法．数值稳定的算法，误差不会放大，甚至会缩小；而数值不稳定的算法会放大误差． 5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab6. 实验过程分别用数组IA[ ]和IB[ ]保存两种算法计算的结果． 7. 结果与分析 运行结果：（或拷屏）8. 附录：程序清单%*************************************************************%* 程序名：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=1:11fprintf('%2d %14.6f %14.6f %14.6f\n',n-1,IA(n),IB(n),I(n));end% n显示为2位整数, 其它显示为14位其中小数点后显示6位的小数《数值计算方法》实验1报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验1 算法设计原则（除数绝对值不能太小） 2. 实验题目将线性方程组增广矩阵利用初等行变换可化为⎪⎪⎭⎫⎝⎛→-⎪⎪⎭⎫ ⎝⎛→-⎪⎪⎭⎫ ⎝⎛''0'0''02221112'12221121112222211121122121121b a b a r r b a b a a r r b a a b a a a a a a由此可解得'/',/'22221111a b x a b x ==．分别解增广矩阵为161011212-⎛⎫ ⎪⎝⎭和162121011-⎛⎫⎪⎝⎭的方程组，验证除数绝对值远小于被除数绝对值的除法会导致结果失真． 3. 实验目的验证数值算法需遵循的若干规则． 4. 基础理论设计数值算法时，应避免除数绝对值远小于被除数绝对值的除法，否则绝对误差会被放大，使结果失真． 5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab6. 实验过程用二维数组A 和B 存放方程组的增广矩阵，利用题目所给初等行变换求解方程组． 7. 结果与分析第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]; % 增广矩阵% 方程组A% 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行加第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))]);《数值计算方法》实验2报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验2 非线性方程的迭代解法(之简单迭代法) 2. 实验题目用简单迭代法求方程010423=-+x x 在区间[1,2]的一个实根，取绝对误差限为410-．3. 实验目的掌握非线性方程的简单迭代法． 4. 基础理论简单迭代法：将方程0)(=x f 改写成等价形式)(x x ϕ=，从初值0x 开始，使用迭代公式)(1k k x x ϕ=+可以得到一个数列，若该数列收敛，则其极限即为原方程的解．取数列中适当的项可作为近似解． 5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程7. 结果与分析8. 附录：程序清单《数值计算方法》实验2报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验2 非线性方程的迭代解法(之Newton 迭代法) 2. 实验题目用Newton 迭代法求方程010423=-+x x 在区间[1,2]的一个实根，取绝对误差限为410-．3. 实验目的掌握求解非线性方程的Newton 迭代法． 4. 基础理论Newton 迭代法：解方程0)(=x f 的Newton 迭代公式为)(')(1k k k k x f x f x x -=+．5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程7. 结果与分析8. 附录：程序清单《数值计算方法》实验2报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验2 非线性方程的迭代解法(之对分区间法) 2. 实验题目用对分区间法求方程310x x --=在区间[1, 1.5]的一个实根，取绝对误差限为410-． 3. 实验目的掌握求解非线性方程的对分区间法． 4. 基础理论对分区间法：取[a ，b ]的中点p ，若f (p ) ≈ 0或b – a < ε，则p 为方程0)(=x f 的近似解；若f (a ) f (p ) < 0，则说明根在区间取[a ，p ]中；否则，根在区间取[p ，b ]中．将新的有根区间记为 [a 1，b 1]，对该区间不断重复上述步骤，即可得到方程的近似根． 5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程用宏定义函数f (x )；为了循环方便，得到的新的有根区间始终用[a ，b ]表示；由于新的有根区间可能仍以a 为左端点，这样会反复使用函数值f (a )，为减少运算次数，将这个函数值保存在一个变量fa 中；同样在判断新的有根区间时用到函数值f (p )，若新的有根区间以p 为左端点，则下一次用到的f (a )实际上就是现在的f (p )，为减少运算次数，将这个函数值保存在一个变量fp 中．算法的伪代码描述：Input ：区间端点a ，b ；精度要求（即误差限）ε；函数f （x ）；最大对分次数N Output ：近似解或失败信息7. 结果与分析8. 附录：程序清单说明: 源程序中带有数字的空行，对应着算法描述中的行号%**********************************************************%* 程序名：Bisection.m *%* 程序功能：使用二分法求解非线性方程. *%**********************************************************f=inline('x^3-x-1'); % 定义函数f(x)a=input('有根区间左端点: a=');b=input('右端点:b=');epsilon=input('误差限:epsilona=');N=input('最大对分次数: N=');1 % 对分次数计数器n置12 % 左端点的函数值给变量fafprintf('\n k p f(p) a(k) f(a(k))'); fprintf(' b(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% while n≤ N 4 % 取区间中点p5% 求p 点函数值给变量fpfprintf('%2d%12.6f%12.6f',n,p,fp); % 输出迭代过程中的中点信息p 和f(p)6 % 如果f(p)=0或b-a 的一半小于误差限εfprintf('\n\n 近似解为:%f\n',p);% 则输出近似根p (7)return;% 并完毕程序 (7)89 % 计数器加110% 若f(a)与f(p)同号11% 则取右半区间为新的求根区间, 即a 取作p 12 % 保存新区间左端点的函数值 13% 否则14 % 左半区间为新的求根区间, 即b 取作p 15fprintf('%12.6f%12.6f%12.6f%12.6f\n',a,fa,b,b-a); %显示新区间端点与左端函数值、区间长度 16fprintf('\n\n 经过%d 次迭代后未达到精度要求.\n',N); % 输出错误信息(行17)《数值计算方法》实验2报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验2 非线性方程的迭代解法(之Aitken-Steffensen 加速法) 2. 实验题目用Aitken-Steffensen 加速法求方程010423=-+x x 在区间[1,2]的一个实根，取绝对误差限为410-．3. 实验目的熟悉求解非线性方程的Aitken-Steffensen 加速法． 4. 基础理论将方程0)(=x f 改写成等价形式)(x x ϕ=，得到从初值0x 开始的迭代公式)(1k k x x ϕ=+后，基于迭代公式)(1k k x x ϕ=+的Aitken-Steffensen 加速法是通过“迭代-再迭代-加速”完成迭代的，具体过程为kk k k k k k k k k k x y z z y x x y z x y +---===+2)(),(),(21ϕϕ． 5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程为了验证Aitken-Steffensen 加速法可以把一些不收敛的迭代加速成迭代收敛，我们使用将方程组变形为31021x x -=，取迭代函数31021)(x x -=ϕ，并利用宏定义出迭代函数．由于不用保存迭代过程，所以用x0表示初值同时也存放前一步迭代的值，y 和z 是迭代过程中产生的y k 和z k ，x 存放新迭代的结果．算法的伪代码描述：Input ：初值x 0；精度要求（即误差限）ε；迭代函数φ（x ）；最大迭代次数N7. 结果与分析8. 附录：程序清单%************************************************************* %* 程 序 名：Aitken_Steffensen.m * %* 程序功能：用Aitken-Steffensen 加速法求方程. * %************************************************************* clc;clear all;phi=inline('0.5 * sqrt( 10 - x^3)'); % 迭代函数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位显示x01 % n 是计数器2 % while n<=Ny= 3 ; % 迭代 z= 3 ; % 再迭代 x= 3 ; % 加速% x0初值与前一步的近似值, y 和z 是中间变量, x 是下一步的近似值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报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验2 非线性方程的迭代解法（之Newton 下山法） 2. 实验题目用Newton 下山法求方程010423=-+x x 在区间[1,2]的一个实根，取绝对误差限为410-．3. 实验目的熟悉非线性方程的Newton 下山法． 4. 基础理论Newton 下山法：Newton 下山法公式为)(')(1k k kk k x f x f x x λ-=+，使|)(||)(|1k k x f x f <+，其中10≤<k λ．5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程定义函数f(x)和df(x)，其中df(x)是f(x)的导函数．每步迭代时先取下山因子为1，尝试迭代，判断尝试结果是否满足下山因子，若满足则作为这步的迭代结果；否则将下山因子减半，然后再尝试．为防止当前的x k 是极小值点，附近不会有满足下述条件的其它点，使尝试陷入死循环，同时计算机中能表示出的浮点数也有下界，因此我们设置了最大尝试次数．当超过最大尝试次数时，不再进行下山尝试．由于反复尝试迭代且要判断下山条件，所以f (x 0)和f ‘(x 0)会反复使用，为避免重复计算浪费运行时间，将这两个值分别保存在变量fx0和dfx0．而尝试产生的节点，判断下山条件时要用到它的函数值，若尝试成功，这个点会作为下一步的初值再使用，所以把该点的函数值也保存在变量fx 中．算法的伪代码描述：Input ：初值x 0；精度要求（即误差限）ε；函数与其导函数f （x ）和f’(x)；最大迭代次数N ；K 下山尝试最大次数Output ：近似解或失败信息7. 结果与分析8. 附录：程序清单%*************************************************************%* 程序名：NewtonDownhill.m *%* 程序功能：用Newton下山法求解非线性方程. *%*************************************************************clc;clear all;f=inline('x^3-x-1'); % 函数f(x)df=inline('3*x^2-1'); % 函数f(x)的导函数x0=input('初值: x0 = ');epsilon=input('误差限: epsilon=');N=input('最大迭代次数: N=');K=input('最大下山尝试次数: K=');1 % 迭代次数计数器2 % 存x0点函数值fprintf('\n\n n x(n) f(x(n))\n'); % 显示表头fprintf('%2d%14.6f%14.6f\n',0,x0,fx0); % 2位整数显示0, 共14位小数6位显示x0和fx03 % while n≤ Ndisp(''); % 换行显示下山尝试过程的表头disp(' 下山因子尝试x(n) 对应f(x(n)) 满足下山条件');disp('');4 % 存x0点导数值, 每次下山尝试不用重新计算ifdfx0==0 % 导数为0不能迭代disp(‘无法进行Newton迭代’);return;endlambda=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)<abs(fx0)) % 判断是否满足下山条件fprintf(' 满足\n');break; % 是, 则退出下山尝试的循环elsefprintf(' 不满足\n');endlambda=lambda/2; % 不是, 则下山因子减半k=k+1; % 计数器加1endif k>Kfprintf('\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\n',x); % (23)return; % 达到, 则显示结果并完毕程序(23) end % (24)% 用x0,fx0存放前一步的近似值和它的函数值, 进行循环迭代25262728fprintf('\n 迭代%d次还不满足误差要求.\n\n',N);《数值计算方法》实验2报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验2 非线性方程的迭代解法（之弦截法） 2. 实验题目用弦截法求方程010423=-+x x 在区间[1,2]的一个实根，取绝对误差限为410-． 3. 实验目的熟悉非线性方程的弦截法． 4. 基础理论将Newton 迭代法中的导数用差商代替，得到弦截法(或叫正割法)公式)()()(111k k k k k k k x f x f x f x x x x --+---=．5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程不保存迭代过程，所以始终以x 0和x 1分别存放x k -1和x k ，而x 存放新产生的迭代值x k +1，这样，下一次迭代时需要把上一步的x 1（即x k ）赋值于x 0（做新的x k -1）．这些点的函数值会重复用到，在迭代公式中也要用到，上一步的x 1作为下一步的x 0也会再一次用它的函数值，为减少重新计算该点函数值的运行时间，将x 1点的函数值保存在变量fx1中．算法的伪代码描述：Input ：初值x 0，x 1；精度要求（即误差限）ε；函数f （x ）；最大迭代次数N7. 结果与分析8. 附录：程序清单%*************************************************************%* 程序名：SecantMethod.m *%* 程序功能：用弦截法求解非线性方程. *%*************************************************************clc;clear all;f=inline('2*x^3-5*x-1'); % 函数f(x)x0=input('第一初值: x0 = ');x1=input('第二初值: x1 = ');epsilon=input('误差限: epsilon=');N=input('最大迭代次数: N=');fprintf('\n n x(n)\n'); % 显示表头fprintf('%2d%14.6f\n', 0, x0); % 占2位显示步数0, 共14位其中小数6位显示x0fprintf('%2d%14.6f\n', 1, x1); % 占2位显示步数1, 共14位其中小数6位显示x11 % 存x0点函数值2 % 存x1点函数值3 % 迭代计数器4 % while n≤ N% 弦截法公式fprintf('%2d%14.6f\n', n, x); %显示迭代过程6 % 达到精度要求否fprintf('\n\n 方程的近似解为: x≈%f\n\n', x);return; % 达到, 则显示结果并完毕程序89 % 原x1做x0为前两步的近似值10 % 现x做x1为一两步的近似值11 % x0点函数值12 % 计算x1点函数值, 为下一次循环13 % 计数器加1 14fprintf('\n 迭代%d 次还不满足误差要求.\n\n',N);《数值计算方法》实验3报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验3 解线性方程组的直接法(之Gauss 消去法) 2. 实验题目用Gauss 消去法求解线性方程组⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--000.3000.2000.1643.5072.1000.2623.4712.3000.1000.3000.2001.0321x x x ． 3. 实验目的掌握解线性方程组的Gauss 消去法． 4. 基础理论Gauss 消去法是通过对增广矩阵的初等行变换，将方程组变成上三角方程组，然后通过回代，从后到前依次求出各未知数．Gauss 消去法的第k 步(1≤k≤n -1)消元：若0≠kk a ，则依次将增广矩阵第k 行的kk ik a a /-倍加到第i 行（k+1≤i≤n），将第k 列对角线下的元素都化成0．5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程7. 结果与分析8. 附录：程序清单《数值计算方法》实验3报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验3 解线性方程组的直接法(之Gauss 列主元消去法) 2. 实验题目用Gauss 列主元消去法求解线性方程组⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--000.3000.2000.1643.5072.1000.2623.4712.3000.1000.3000.2001.0321x x x ． 3. 实验目的掌握解线性方程组的Gauss 列主元消去法． 4. 基础理论Gauss 列主元消去法也是通过对增广矩阵的初等行变换，将方程组变成上三角方程组，然后通过回代，从后到前依次求出各未知数．Gauss 列主元消去法的第k 步(1≤k≤n -1)消元：先在nk k k kk a a a ,,,,1 +中找绝对值最大的，将它所在的行与第k 行交换，然后将第k 行的kk ik a a /-倍加到第i 行（k+1≤i≤n），将第k 列对角线下的元素都化成0． 5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程7. 结果与分析8. 附录：程序清单《数值计算方法》实验3报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验3 解线性方程组的直接法(之Doolittle 分解) 2. 实验题目对矩阵A 进行Doolittle 分解，其中⎪⎪⎪⎪⎪⎭⎫⎝⎛----=3101141101421126A ．3. 实验目的掌握矩阵的Doolittle 分解． 4. 基础理论矩阵的Doolittle 分解是指将矩阵n n ij a A ⨯=)(可以分解为一个单位下三角矩阵和一个上三角矩阵的乘积．若设⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=nn n n n n n n u u u u u u u u u u U l l ll l l L000000,1010010001333223221131211321323121则可依如下顺序公式计算⎪⎪⎩⎪⎪⎨⎧++=-=+=-=∑∑-=-=1111,,2,1,/)(,,1,,k t kk tk it ik ik k r rj kr kj kj nk k i u u l a l nk k j u l a u其中k = 1，2，…，n ．5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程（1）按计算公式依次计算一行u 同时计算一列l ；（2）因为计算完u ij （或l ij ）后，a ij 就不再使用，为节省存储空间，将计算的u ij （和l ij ）仍存放在矩阵A 中的相应位置；（3）使用L 矩阵和U 矩阵时需要根据元素所在位置取固定值或A 中相应位置的值．L 对角线上的元素为1，上三角部分为0，下三角部分为A 中对应的元素；U 的下三角部分为0，上三角部分为A 中对应的元素．算法的伪代码描述： Input ：阶数n ；矩阵A7. 结果与分析8. 附录：程序清单%****************************************************% 程序名: Doolittle.m *% 程序功能: 矩阵LU分解中的Doolittle分解. *%****************************************************clc;clear all;n=4; % 矩阵阶数A=[6 2 1 -1;2 4 1 0; 1 1 4 -1; -1 0 -1 3]disp('A=');disp(A);% LU分解(Doolittle分解)for k=1:n% 计算矩阵U的元素u_{kj}% (可参照下面l_{ik}的公式填写)% 计算矩阵L的元素l_{ik}% L 在A 下三角, U 在上三角(对角线为1) enddisp('分解结果：'); disp('L='); for i=1:n for j=1:nif i>j % 在下三角部分, 则取A 对于的元素显示 fprintf(' %8.4f',A(i,j));elseif i==j % 在对角线上, 则显示1 fprintf(' %8d',1);else % 在上三角部分, 则显示0 fprintf(' %8d',0); end endfprintf('\n'); % 换行 enddisp('U='); for i=1:n for j=1:nif i<=j % 在上三角部分或对角线上, 则取A 对于的元素显示 fprintf(' %8.4f',A(i,j));else % 在下三角部分, 则显示0 fprintf(' %8d',0); end endfprintf('\n'); % 换行 end《数值计算方法》实验3报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验3 解线性方程组的直接法(之LU 分解法) 2. 实验题目用LU 分解(Doolittle 分解)法求解线性方程组⎪⎩⎪⎨⎧=++=++=++104615631552162321321321x x x x x x x x x 3. 实验目的熟悉解线性方程组LU 分解法．4. 基础理论若将矩阵A 进行了Doolittle 分解，A = LU ，则解方程组b x A=可以分解求解两个三角方程组b y L=和y x U =．它们都可直接代入求解，其中b y L=的代入公式为∑-==-=11,,2,1,k j j kj k k n k y l b y而y x U=的代入公式为∑+=-=-=nk j kk j kjk k n n k u x uy x 11,,1,,/)( ．5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程（1）Doolittle 分解过程依次计算一行u 同时计算一列l 完成，并将计算的u ij （和l ij ）仍存放在矩阵A 中的相应位置；（2）求解方程组的代入公式中用到的u ij 和l ij 都直接在A 的相应位置取值即可． 算法的伪代码描述：Input ：阶数n ；矩阵A ；常数项向量b7. 结果与分析8. 附录：程序清单%**************************************************** % 程序名: LinearSystemByLU.m *% 程序功能: 利用LU分解(Doolittle分解)解方程组. *%****************************************************clc;clear all;n=3; % 矩阵阶数A=[1 2 6; 2 5 15; 6 15 46];b=[1;3;10];% LU分解(Doolittle分解)for k=1:n% 计算矩阵U的元素u_{kj}% (可参照下面l_{ik}的公式填写)% 计算矩阵L的元素l_{ik}% L在A下三角, U在上三角(对角线为1) endfor k=1:n % 用代入法求解下三角方程组Ly=by(k)=b(k);3 %∑-==-=11,,2,1,kjj kjk knkylby33enddisp('方程组Ly=b的解：y=');disp(y');for k=n:-1:1 % 回代求解上三角方程组Ux=y x(k)=y(k);6 %∑+=-=-=nkjj kjk knnkxuyx11,,1,,666 enddisp('原方程组的解：x='); disp(x');《数值计算方法》实验3报告班级： 20##级####x 班 学号： 20##2409#### ： ##X成绩：1. 实验名称实验3 解线性方程组的直接法(之Cholesky 分解) 2. 实验题目对矩阵A 进行Cholesky 分解，其中⎪⎪⎪⎪⎪⎭⎫⎝⎛----=3101141101421126A ． 3. 实验目的理解矩阵的Cholesky 分解． 4. 基础理论矩阵的Cholesky 分解是指将矩阵n n ij a A ⨯=)(可以分解为一个下三角矩阵L 和L 转置的乘积，即A =LL T，其中L 各元素可依如下顺序公式计算⎪⎪⎩⎪⎪⎨⎧++=-=-=∑∑-=-=11112,,2,1,/)(k t kktk it ik ik k r kr kk kk nk k i l l l a l l a l其中k = 1，2，…，n ．5. 实验环境操作系统：Windows xp ； 程序设计语言：VC++ 6. 实验过程（1）按计算公式依次先计算一列对角线上的元素l kk ，再计算这列其他元素l ik ，且对称位置的元素也取同一个值；（2）因为计算完l ij 后，a ij 就不再使用，为节省存储空间，将计算的l ij 仍存放在矩阵A 中的相应位置；（3）使用L 矩阵时需要根据元素所在位置取固定值或A 中相应位置的值．L 上三角部分为0，对角线和下三角部分为A 中对应的元素．算法的伪代码描述：Input ：阶数n ；矩阵AOutput ：矩阵L (合并存储在数组A 中)行号 伪代码注释1 for k ← 1 to n2∑-=-=112k r krkk kk l a l3 for i ← k to n4 ∑-=-=11/)(k t kk tk it ik ik l l l a l计算结果存放在a ij5 endfor6 endfor7return L输出L7. 结果与分析8. 附录：程序清单%************************************************************* %* 程 序 名：Cholesky.m * %* 程序功能：对称正定矩阵的Cholesky 分解. * %*************************************************************n=4; % 矩阵阶数 A=[6,2,1,-1; 2,4,1,0; 1,1,4,-1; -1,0,-1,3];disp('A ='); for i=1:n for j=1:nfprintf('%10.4f',A(i,j)); % 共占14位endfprintf('\n');% 一行完毕换行end% Cholesky 分解 for k=1:n % 计算对角线上的l _{kk}% 计算其他的l _{ik} % 和l _{ki}end % L 在A 下三角, L^T 在上三角disp('分解结果：'); disp('L='); for i=1:n for j=1:n if i>=j % 在下三角部分或对角线上, 则取A 对于的元素显示fprintf('%10.4f',A(i,j));else % 在上三角部分, 则显示0 fprintf('%10d',0); end endfprintf('\n'); % 换行 end《数值计算方法》实验3报告班级： 20##级####x 班 学号： 20##2409#### ： ##X成绩：1. 实验名称实验3 解线性方程组的直接法(之改进的Cholesky 分解) 2. 实验题目对矩阵A 进行改进的Cholesky 分解，其中⎪⎪⎪⎪⎪⎭⎫⎝⎛----=3101141101421126A ．3. 实验目的理解矩阵改进的Cholesky 分解． 4. 基础理论矩阵的改进的Cholesky 分解是指将矩阵n n ij a A ⨯=)(可以分解为一个单位下三角矩阵L 和对角矩阵D 与L 转置的乘积，即A =LDL T，其中L 和D 各元素可依如下顺序公式计算⎪⎪⎩⎪⎪⎨⎧++=-=-=∑∑-=-=11112,,2,1,/)(k t k kt it t ik ik k r kr r kk k nk k i d l l d a l l d a d其中k = 1，2，…，n ．5. 实验环境操作系统：Windows xp ； 程序设计语言：VC++ 6. 实验过程（1）按计算公式依次先计算D 的一个元素d k ，再计算L 中这列的元素l ik ，且对称位置的元素也取同一个值；（2）因为计算完d k 和l ij 后，a kk 或a ij 就不再使用，为节省存储空间，将计算的a kk 或l ij 仍存放在矩阵A 中的相应位置；（3）使用L 矩阵时需要根据元素所在位置取固定值或A 中相应位置的值．L 对角线和上三角部分为0，下三角部分为A 中对应的元素；D 对角线为A 中对应的元素，其余都是0．算法的伪代码描述： Input ：阶数n ；矩阵AOutput ：矩阵L (合并存储在数组A 中)7. 结果与分析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 =');for i=1:nfor j=1:nfprintf('%10.4f',A(i,j)); % 共占14位endfprintf('\n'); % 一行完毕换行end% Cholesky分解for k=1:n% 计算D对角线上的u_{kk}% 计算L的元素l_{ik}% 和L转置的元素l_{ki} end % L在A下三角, D在对角线disp('分解结果：');disp('L=');for i=1:nfor j=1:nif i>j % 在下三角部分, 则取A对于的元素显示fprintf('%10.4f',A(i,j));elseif i==j % 在对角线上, 则显示1fprintf('%10d',1);else % 在上三角部分, 则显示0fprintf('%10d',0);endendfprintf('\n'); % 换行enddisp('D='); for i=1:n for j=1:n if i==j % 在对角线上, 则取A 对于的元素显示fprintf('%10.4f',A(i,j));else % 其余显示0fprintf('%10d',0); end endfprintf('\n'); % 换行 end《数值计算方法》实验3报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验3 解线性方程组的直接法(之追赶法) 2. 实验题目用追赶法求解线性方程组⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-----101053001210023100124321x x x x 3. 实验目的熟悉解线性方程组的追赶法． 4. 基础理论对于系数矩阵为三对角矩阵的方程组，其Crout 分解可分解为⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=------11111211122111122211n n nn n n nn n n t t t s a s a s a s b a c b a c b a c b A这样，解方程组可以由如下2步完成：“追”：,,,3,2,/)(,,/,/,1111111111n i s y a f y t a b s s c t s f y b s i i i i i i i i i i i i =-=-====-----其中：Tn f f ),,(1 为方程组的常数项，n t 没用；“赶”：.1,,2,1,,1 --=-==+n n i x t y x y x i i i i n n5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程在“追”的过程中，向量s 和y 都有n 个元素，t 只有n -1个元素，又1s 和1y 的计算公式与其它i s 和i y 不同，所以先单独计算1s 和1y ，然后在一个n -1次循环中，求其它i s 和i y 以与i t ．由于在“追”的过程中，i b ，i c 和i f 在分别计算完对应的i s ，i t 和i y 后就不再使用，所以借用数组b ，c 和f 存储向量s ，t 和y ；同样在“赶”的过程中，i y 在计算完对应的i x 后就不再使用，所以再一次借用数组f 存储向量x ．追赶法算法的伪代码描述：Input ：阶数n ；三对角矩阵的三条对角线向量a ，b ，c ，常数项向量f Output ：方程组的解x改进的追赶法算法的伪代码描述：Input ：阶数n ；三对角矩阵的三条对角线向量a ，b ，c ，常数项向量f Output ：方程组的解x7. 结果与分析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_1 for 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% 再求x_i(i=n-1,n-2, (1)endx=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(n); % x_n等于y_n仍在f_n中for k=n-1:-1:1% 再求x_i(i=n-1,n-2, (1)endx=f' % 输出解向量《数值计算方法》实验4报告班级：20##级####x班学号：20##2409####：##X 成绩：1. 实验名称实验4 解线性方程组的迭代法(之Jacobi迭代)2. 实验题目用Jacobi迭代法求解线性方程组1231231232251223x x x x x x x x x +-=⎧⎪++=⎪⎨++=⎪⎪⎩任取3. 实验目的掌握解线性方程组的Jacobi 迭代法． 4. 基础理论将第i （n i ≤≤1）个方程i n in i i b x a x a x a =+++ 2211移项后得到等价方程ii n in i i i i i i i i i a x a x a x a x a b x /)(11,11,11------=++--便可构造出Jacobi 迭代公式,1,0,/)()()(11,)(11,)(11)1(=------=++--+k a x a x a x a x a b x ii k n in k i i i k i i i k i i k i ． 5. 实验环境操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程7. 结果与分析8. 附录：程序清单《数值计算方法》实验4报告班级： 20##级####x 班 学号： 20##2409#### ： ##X 成绩：1. 实验名称实验4 解线性方程组的迭代法(之Gauss-Seidel 迭代) 2. 实验题目用Gauss-Seidel 迭代法求解线性方程组。
线性方程组的数值解法-安振华-2012011837
实验5：线性方程组的数值解法化学工程系分2 安振华2012011837【实验目的】1、掌握线性方程组的常用数值解法，包括高斯消去法、LU分解法以及校正法。
2、体验数值计算的时间复杂度和计算规模的关系。
3、加深对数值计算误差的理解。
4、学习使用迭代法等算法，求解非线性方程。
5、学习如何使用MATLAB解非线性方程组和方程组。
【实验容】【实验五：习题9】种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应保持不变，种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n，当年年龄k的种群数量记作x k，繁殖率记作b k（每个雌性个体在1年繁殖的数量），自然存活率记作s k（s k=1-d k，d k为1年的死亡率），收获量记作h k，则来年年龄k的种群数量k x应为：111,(1,2,,1)n k k k k k k k x b x x s x h k n +===-=⋅⋅⋅-∑要求各个年龄的种群数量每年维持不变就是要使(1,2,,)k k x x k n ==⋅⋅⋅（1） 若b k ，s k 已知，给定收获量h k ，建立求各年龄的稳定种群数量x k 的模型（用矩阵向量表示）（2） 设n=5，b 1=b 2=b 5=0，b 3=5，b 4=3，s 1=s 4=0.4，s 2=s 3=0.6，如果要求h 1~h 5为500,400,200,100,100,求x 1~x 5 （3） 要使h 1~h 5均为500，如何达到？ 【分析】为方便起见以下种群数量均指其中的雌性。
我们并且有以下的假设：(1)雌性个体的繁殖率和存活率在特定的时间是不变的。
(2)人工饲养的种群在质量和数量上是不受外界环境和资源的限制的。
(3)模型中不考虑人为的或是自然的灾害所造成的种群数量、繁殖率和存活率的变动。
数值计算(分析)实
数值计算(分析)实验报告2南昌航空大学数学与信息科学学院实验报告课程名称：《数值计算方法》实验名称：曲线拟合实验类型：验证性■综合性□设计性□实验室名称：数学实验室班级学号： 09072113学生姓名：邢宪平任课教师（教师签名）：成绩：一、实验目的实验目的：实验目的：了解函数逼近与曲线拟合的基本原理，并且运用MATLAB 软件进行实践操作。
二、实验原理、程序框图、程序代码等 实验题目：题目1：试分别用抛物线2y a bx cx =++和指数曲线bxy ae =拟合下列数据并比较两个拟合函数的优劣。
题目2：已知实验数据如下：试用形如2y a bx =+的抛物线进行最小二乘拟合。
实验原理：1、逼近方式 假设()[,]f x C a b ∈，2{1,,,...,}n nHspan x x x =，()nnP x H ∈，称(,)|||||()()|max n n n a x bf P F P f x P x ≤≤=-=-V 为()f x 与()|nP x 在[,]a b 上的偏差。
若存在*()nnP x H ∈，使得**(,)|||||()()|max inf n nn nn P H a x bf P f Pf x P x ∞∈≤≤=-=-V 则称*()nP x 是()f x 在[,]a b 上的最佳一致逼近多项式。
假设()[,]f x C a b ∈及[,]C a b 的一个子集01{(),(),,...()}nspan x x x ϕ=ϕϕϕ，若存在*()S x ϕ∈，使*22222()()||()()||||()()||()[()()]min min bS x S x af x S x f x S x x f x S x dxϕϕρ∈∈-=-=-⎰则称*()S x 是()f x 在子集[,]C a b ϕ⊂中的最佳平方逼近数。
2、曲线拟合上述函数的最佳平方逼近法中，若()f x 是以一组离散点集的形式给出的，即给出了函数()f x 在一些离散点上的值{(,),0,1,...,}iix y i m =，则该方法就是所说的曲线拟合。
山东大学数值计算实验报告3
2.1.2 直接解方程组
fprintf('%s\n','n=2 时，解法 2');
x12 = H1\b1
运行结果如下：
2.2 当 n=5 时
%生成 n=5 时希尔伯特矩阵 H2=zeros(5,5); n2=5; for i = 1:5 for j = 1:5 H2(i,j) = 1/(i+j-1); end end fprintf('%s\n','n=5 时，解法 1：'); L2=c_1(H2,n2)%楚列斯基函数得到 L %得到 b t2=ones(5,1); b2=H2*t2; %利用之前实验 2 所写的前代，后代函数得到答案 front2 = a_3_a( L2,b2 ); x21 = a_3_b(L2',front2)
%当 n=10 时 fprintf('%s\n','当 n=10 时，该矩阵的 1 范数及 1 范数下的条件数为：'); fan31 = norm(H3,1) cond31 = cond(H3,1) fprintf('%s\n','当 n=10 时，该矩阵的 2 范数及 2 范数下的条件数为：'); fan32 = norm(H3,2) cond32 = cond(H3) fprintf('%s\n','当 n=10 时， 该矩阵的无穷范数及无穷范数下的条件数为： '); fan33 = norm(H3,inf) cond33 = cond(H3,inf)
运行结果如下：
fprintf('%s\n','n=10 时，解法 2'); x32 = H3\b3
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值实验课试题本次数值实验课结课作业，请按题目要求内容写一篇文章。
按题目要求人数自由组合，每组所选题目不得相同(有特别注明的题目除外)。
试题如下：1）解线性方程组的Gauss 消去法和列主元Gauss 消去法(2人)/*张思珍，巩艳华*/用C 语言将不选主元和列主元Gauss 消去法编写成通用的子程序，然后用你编写的程序求解下列84阶的方程组⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛141515151576816816816816816848382321 x x x x x x参考书目：1.《计算机数值方法》，施吉林、刘淑珍、陈桂芝编2.《数值线性代数》，徐树方、高立、张平文编3.《数值分析简明教程》，王能超编2）解线性方程组的平方根法(4人)/*朱春成、黄锐奇、张重威、章杰*/ 用C 语言将平方根法和改进的平方根法编写成通用的子程序，然后用你编写的程序求解对称正定方程组b Ax =，其中(1)b 随机的选取，系数矩阵为100阶矩阵⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛1011101110111011101110 ； (2)系数矩阵为40阶的Hilbert 矩阵，即系数矩阵A 的第i 行第j 列元素为11-+=j i a ij ，向量b 的第i 个分量为∑=-+=n j i j i b 111. 参考书目：1.《计算机数值方法》，施吉林、刘淑珍、陈桂芝编2.《数值线性代数》，徐树方、高立、张平文编3.《数值分析简明教程》，王能超编3）三对角线方程组的追赶法(3人)/*黄佳礼、唐伟、韦锡倍*/用C 语言将三对角线方程组的追赶法法编写成通用的子程序，然后用你编写的程序求解如下84阶三对角线方程组⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛141515151576816816816816816848382321 x x x x x x参考书目：1.《计算机数值方法》，施吉林、刘淑珍、陈桂芝编2.《数值分析简明教程》，王能超编4）线性方程组的Jacobi 迭代法(3人)/*周桂宇、杨飞、李文军*/用C 语言将Jacobi 迭代法编写成独立的子程序，并用此求解下列方程组，精确到小数点后5位⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛-1490122111221321x x x参考书目：1.《计算机数值方法》，施吉林、刘淑珍、陈桂芝编2.《数值线性代数》，徐树方、高立、张平文编3.《数值分析简明教程》，王能超编5）线性方程组的Gauss-Seidel 迭代法(3人)/*张玉超、范守平、周红春*/ 用C 语言将Gauss-Seidel 迭代法编写成独立的子程序，并用此求解下列方程组，精确到小数点后5位⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--397211111112321x x x参考书目：1.《计算机数值方法》，施吉林、刘淑珍、陈桂芝编2.《数值线性代数》，徐树方、高立、张平文编3.《数值分析简明教程》，王能超编6）解线性方程组的最速下降法法(2人)/*赵育辉、阿热孜古丽*/用C 语言将最速下降法编写成通用的子程序，然后用你编写的程序求解对称正定方程组b Ax =，其中b 随机的选取，系数矩阵为100阶矩阵⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛1011101110111011101110 参考书目：1.《数值线性代数》，徐树方、高立、张平文编2.《最优化方法及其应用》，郭科、陈聆、魏友华编7）解线性方程组的共轭梯度法(3人)/*刘森林、钟荣生、芦佩*/用C 语言将共轭梯度法编写成通用的子程序，然后用你编写的程序求解对称正定方程组b Ax =，其中系数矩阵为40阶的Hilbert 矩阵，即系数矩阵A 的第i 行第j 列元素为参考书目：1.《数值线性代数》，徐树方、高立、张平文编2.《最优化方法及其应用》，郭科、陈聆、魏友华编8）Newton 法求多元二次方程的最优值(3人)/*李馨蕾、杨宏宇、李敏*/ 用C 语言将求解多元二次方程最优值的Newton 法编写成通用的子程序，并用此求解2122212141060)(x x x x x x X f -++--=的极小值，初始点为T X ]00[0=，精度为00001.0=ε.参考书目：《最优化方法及其应用》，郭科、陈聆、魏友华编9）共轭梯度法法求多元二次方程的最优值(3人)/*张南、佟雪、杨坤*/ 用C 语言将求解多元二次方程最优值的共轭梯度法编写成通用的子程序，并用此求解22214)(x x X f += 的极小值，初始点为T X ]11[0=，精度为00001.0=ε.参考书目：《最优化方法及其应用》，郭科、陈聆、魏友华编10）求多元二次方程的最优值的变尺度法(DFP 法)(4人)/*史建国、黄嘉莹、方芳、李念超*/用你熟悉的计算机语言将DFP 变尺度法编写成通用的子程序，并用此求解2221)6()5(4)(-+-=x x X f的极小值点，初始点为T X ]98[0=，精度为00001.0=ε.参考书目：《最优化方法及其应用》，郭科、陈聆、魏友华编11）求多元二次方程的最优值的变尺度法(BFGS 法)(4人)/*袁雪华、孙婷婷、郭良、陈乾*/用你熟悉的计算机语言将BFGS 变尺度法编写成通用的子程序，并用此求解2221)6()5(4)(-+-=x x X f的极小值点，初始点为T X ]98[0=，精度为00001.0=ε.参考书目：《最优化方法及其应用》，郭科、陈聆、魏友华编12）QR 法求解最小二乘问题(4人)/*付为政、董泽尧、黄自鹏、武继飞*/ 用C 语言编写利用QR 分解求解线性最小二乘问题的通用子程序，并用此求解一个二次多项式c bt at y ++=2是其在最小二范数意义下拟合下列数据参考书目：《数值线性代数》，徐树方、高立、张平文编13）矩阵求逆(4人)/*陈巧、汪恒、陈朝、何义连*/利用你已有的理论知识，使用C 语言编写一个求解给定矩阵的逆的通用子程序，并用此求解矩阵⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=2271.02168.12071.01968.01871.01768.01675.11582.01490.01397.01254.01161.12671.12568.02471.02368.0A 的逆。
参考书目： 1.《计算机数值方法》，施吉林、刘淑珍、陈桂芝编2.《数值线性代数》，徐树方、高立、张平文编3.《数值分析简明教程》，王能超编14）一元非线性方程求根的迭代法(4人)/*周孝金、雍佳飞、史旭吉、罗职权*/用你熟悉的计算机语言将求解一元非线性方程的牛顿法、弦截法和快速弦截法编写成独立的子程序；并分别用此求解下列方程的全部根020102)(23=-++=x x x x f要求精度为0.000001.参考书目：1.《数值分析简明教程》，王能超编2.《计算机数值方法》，施吉林、刘淑珍、陈桂芝编15）无约束最优化问题的一维搜索法(4人)/*覃伟、班光德、吴本远、韦祥胜*/用你熟悉的计算机语言将无约束最优化问题的一维搜索法(包括对分法、黄金分割法、Newton 切线法、抛物线插值法)编写成独立的子程序；并分别用此求3728)(23+--=x x x x f在区间[]]2,0[,=b a 的最小值点，要求精度为0.0001.参考书目：《最优化方法及其应用》，郭科、陈聆、魏友华编16）常微分方程数值解的Adams 法(3人)/*莫尚威、杨令宗、高朝家*/ 用C 语言将常微分方程数值解的Adams 法)编写成独立的子程序；并分别用此求⎪⎩⎪⎨⎧=-=1)0(2'y y x y y 的数值解，取步长为1.0=h .参考书目：《数值分析简明教程》，王能超编17）高震荡函数的数值积分(3人)/*李杨鹏、刁洋、田万忠*/请用你熟悉的计算机语言编写适合高震荡函数的数值积分的独立子程序，并用此求定积分(1)⎰=π5.10115cos xdx I ,(精确值为0666667.0151≈)， (2)⎰=π20230sin cos xdx x x I ,(精确值为209672222.089960-≈-π)，(3)⎰=π20330cos cos xdx x x I ,(精确值为0)； 要求精度为0000001.0=ε.注：此题允许两组同学同时选作。
18）复数运算(3人)/*肖有忠、马成虎、马金云*/ 用C 语言实现复数的加法、减法、乘法、除法、乘幂、n 次方根运算，以及复数指数，复系数多项式的乘积，复矩阵相加、相减、相乘等运算。
19）矩阵的LU 分解(3人)/*田小兵、朱雪梅、祖丽合马*/用你熟悉的计算机语言将矩阵的LU 分解写成通用的子程序，并用此求矩阵 ⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-11242142612332442 的LU 分解；要求输出矩阵L 、U.参考书目： 1.《数值线性代数》，徐树方、高立、张平文编2.《计算机数值方法》，施吉林、刘淑珍、陈桂芝编20）矩阵的QR 分解(4人)用C 语言将矩阵的QR 分解写成通用的子程序，并用此求矩阵⎪⎪⎪⎪⎪⎭⎫ ⎝⎛---121011012111的QR 分解；要求输出矩阵Q 、R.注：一共20题，除第17题可以由两组同学选作外，其他每题只能由一组同学选作。
无特别要求的题目，在编写程序时可选择计算机语言为C 语言、Matlab 、C++、java.文章按格式书写：题目使用三号黑体字、一级标题使用四号黑体字，并居中；二级、三级标题用小四号黑体字，左端对齐（不居中）。
文中其他汉字一律采用小四号宋体字，行距采用1.25倍行距。
文章中至少包括以下六部分内容内容1.算法设计或算法分析2.算法实现(步骤)3.源程序代码4.运算结果5.误差分析6.总结。