浙江大学计算方法大作业第八题
数值分析大作业-三、四、五、六、七
大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) 目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20) 4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28) 数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10） 1.2牛顿法 程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10） 大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =-= ?-???解： Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; 计算方法上机报告 姓名： 学号： 班级： 上课班级： 说明： 本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。 1. 对以下和式计算： ∑ ∞ ? ?? ??+-+-+-+=0681581482184161n n n n S n ，要求： ① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算； （1） 算法思想 1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 1421114 16818485861681 n n n a n n n n n ε??= ---<< ?+++++??； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数) （2）算法结构 1. ;0=s ?? ? ??+-+-+-+= 681581482184161n n n n t n ； 2. for 0,1,2,,n i =??? if 10m t -≤ end; 3. for ,1,2,,0n i i i =--??? ;s s t =+ （3）Matlab源程序 clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0; for n=0:50 t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); if t<=10^(-m) %判断通项与精度的关系break; end end; fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值 for i=n-1:-1:0 t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s=s+t; %求和运算 end s=vpa(s,m) %控制s的精度 （4）结果与分析 当保留11位有效数字时，需要将n值加到n=7， s =3.1415926536； 当保留30位有效数字时，需要将n值加到n=22, s =3.14159265358979323846264338328。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。 参考答案 第一章 1 *1x =1.7； * 2x =1.73； *3x =1.732 。 2． 3． （1） ≤++)(* 3*2*1x x x e r 0.00050； （注意：应该用相对误差的定义去求） （2） ≤)(*3*2*1x x x e r 0.50517； （3） ≤)/(*4*2x x e r 0.50002。 4．设6有n 位有效数字，由6≈2.4494……，知6的第一位有效数字1a =2。 令3)1()1(1* 102 1 102211021)(-----?≤??=?= n n r a x ε 可求得满足上述不等式的最小正整数n =4，即至少取四位有效数字，故满足精度要求可取6≈2.449。 5. 答：（1）*x (0>x )的相对误差约是* x 的相对误差的1/2倍； （2）n x )(* 的相对误差约是* x 的相对误差的n 倍。 6． 根据******************** sin 21)(cos 21sin 21)(sin 21sin 21)(sin 21)(c b a c e c b a c b a b e c a c b a a e c b S e r ++≤ =* *****) ()()(tgc c e b b e a a e ++ 注意当20* π < 7．设20= y ，41.1*0 =y ，δ=?≤--2* 00102 1y y 由 δ1* 001*111010--≤-=-y y y y ， δ2*111*221010--≤-=-y y y y M δ10*991*10101010--≤-=-y y y y 即当0y 有初始误差δ时，10y 的绝对误差的绝对值将减小10 10-倍。而110 10 <<-δ，故计算过程稳定。 8. 变形后的表达式为： （1）)1ln(2--x x =)1ln(2-+-x x （2）arctgx x arctg -+)1(=) 1(11 ++x x arctg （3） 1ln )1ln()1(ln 1 --++=? +N N N N dx x N N =ΛΛ+-+- +3 2413121)1ln(N N N N 1ln )11ln()1(-++ +=N N N N =1)1ln()1 1ln(-+++N N N （4）x x sin cos 1-=x x cos 1sin +=2x tg 计算方法上机报告 姓名: 学号: 班级: 目录 题目一----------------------------------------------------------------- 4 - 1.1题目内容-------------------------------------------------------- 4 - 1.2算法思想-------------------------------------------------------- 4 - 1.3Matlab 源程序----------------------------------------------------- 5 - 1.4计算结果及总结------------------------------------------------- 5 - 题目二----------------------------------------------------------------- 7 - 2.1题目内容-------------------------------------------------------- 7 - 2.2算法思想-------------------------------------------------------- 7 - 2.3 Matlab 源程序---------------------------------------------------- 8 - 2.4计算结果及总结------------------------------------------------- 9 - 题目三--------------------------------------------------------------- -11- 3.1题目内容----------------------------------------------------------- 11 - 3.2算法思想----------------------------------------------------------- 11 - 3.3Matlab 源程序--------------------------------------------------- -13 - 3.4计算结果及总结----------------------------------------------------- 14 - 题目四--------------------------------------------------------------- -15 - 4.1题目内容----------------------------------------------------------- 15 - 4.2算法思想----------------------------------------------------------- 15 - 4.3Matlab 源程序--------------------------------------------------- -15 - 4.4计算结果及总结----------------------------------------------------- 16 - 题目五--------------------------------------------------------------- -18 - -18 - 5.1题目内容 5.2算法思想----------------------------------------------------------- 18 - 5.3 Matlab 源程序--------------------------------------------------- -18 - 2020年奥鹏吉大网络教育《计算方法》大作业解答 (说明：前面是题目，后面几页是答案完整解答部分，注意的顺序。) 一、解线性方程 用矩阵的LU分解算法求解线性方程组 用矩阵的Doolittle分解算法求解线性方程组 用矩阵的Doolittle分解算法求解线性方程组 用高斯消去法求解线性方程组 用高斯消去法求解线性方程组 用主元素消元法求解线性方程组 用高斯消去法求解线性方程组 利用Doolittle分解法解方程组Ax=b，即解方程组 1、用矩阵的LU分解算法求解线性方程组 X1+2X2+3X3 = 0 2X1+2X2+8X3 = -4 -3X1-10X2-2X3 = -11 2、用矩阵的Doolittle分解算法求解线性方程组 X1+2X2+3X3 = 1 2X1– X2+9X3 = 0 -3X1+ 4X2+9X3 = 1 3、用矩阵的Doolittle分解算法求解线性方程组 2X1+X2+X3 = 4 6X1+4X2+5X3 =15 4X1+3X2+6X3 = 13 4、用高斯消去法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 5、用无回代过程消元法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 6、用主元素消元法求解线性方程组 2X 1- X 2+3X 3 = 2 4X 1+2X 2+5X 3 = 4 -3X 1+4X 2-3X 3 = -3 7、用高斯消去法求解线性方程组 123123123234 4272266 x x x x x x x x x -+=++=-++= 8、利用Doolittle 分解法解方程组Ax=b ，即解方程组 12341231521917334319174262113x x x x -? ????? ???? ??-??????=? ? ????--?????? --???? ?? 数值分析报大作业 班级：铁道2班 专业：道路与铁道工程 姓名：蔡敦锦 学号：13011260 一、序言 该数值分析大作业是通过C语言程序编程在Microsoft Visual C++ 6.0编程软件上运行实现的。本来是打算用Matlab软间来计算非线性方程的根的。学习Matlab也差不多有一个多月了，感觉自己编程做题应该没什么问题了；但是当自己真心的去编程、运行时才发现有很多错误，花了一天时间修改、调试程序都没能得到自己满意的结果。所以，我选择了自己比较熟悉的C程序语言来编程解决非线性的求值问题，由于本作业是为了比较几种方法求值问题的收敛速度和精度的差异，选择了一个相对常见的非线性函数来反映其差异，程序运行所得结果我个人比较满意。编写C语言，感觉比较上手，程序出现问题也能比较熟练的解决。最终就决定上交一份C程序语言编程的求值程序了！ 二、选题 本作业的目的是为了加深对非线性方程求根方法的二分法、简单迭代法、、牛顿迭代法弦截法等的构造过程的理解；能将各种方法的算法描述正确并且能够改编为程序并在计算机上实现程序的正确合理的运行，能得到自己满意的结果，并且能调试修改程序中可能出现的问题和程序功能的增减修改。本次程序是为了比较各种方法在求解同一非线性方程根时，在收敛情况上的差异。 为了达到上面的条件我选择自己比较熟悉的语言—C语言来编程，所选题目为计算方程f(x)=x3-2x-5=0在区间[2,3]内其最后两近似值的差的绝对值小于等于5 ?的根的几种方法的比较。 110- 本文将二分法、牛顿法、简单迭代法、弦截法及加速收敛法这五种方法在同一个程序中以函数调用的方式来实现，比较简洁明了，所得结果能很好的比较，便于分析；发现问题和得出结论。 《计算方法B》上机 实验报告 学院：机械工程学院 班级: 姓名： 学号： 2015年12月22日 1 1.计算以下和式： S = ∑ 8n + 1 - 8n + 4 - 8n + 5 - 8n + 6 ? ，要求： 4 2 1 1 ∞ n =0 16 n ? ? ? ? (1)若保留 11 个有效数字，给出计算结果，并评价计算的算法； (2)若要保留 30 个有效数字，则又将如何进行计算。 实现思想： 以上问题出现了近似数相减的问题，为了减小误差，可分别求得减数之和 以及被减数之和，最后将两者相减。另外，减数与被减数求和均为同号计算， 按照绝对值递增顺序相加可减小舍入误差。此题中对有效数字有要求，因而计 算时首先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。 源程序： m=input('输入有效数字个数m='); s0=1;s1=0;s2=0;n=0; %判断迭代次数 while s0>=0.5*10^-(m-1) s0=4/(16^n*(8*n+1))-2/(16^n*(8*n+4))-1/(16^n*(8*n+5))- 1/(16^n*(8*n+6)); n=n+1; end %分别求解各项并求和 for k=n-1:-1:0 a1=4/(16^k*(8*k+1)); a2=2/(16^k*(8*k+4)); a3=1/(16^k*(8*k+5)); a4=1/(16^k*(8*k+6)); s1=a1+s1; s2=a4+a3+a2+s2; end S=vpa(s1-s2,m) 实验结果：11位有效数字计算结果如图1所示；30为有效数字计算结果如图2所示。 图1.11位有效数字计算结果图2.30为有效数字计算结果 计算方法实验报告 班级： 学号： 姓名： 成绩： 1 舍入误差及稳定性 一、实验目的 （1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； （2）通过上机计算，了解舍入误差所引起的数值不稳定性 二、实验内容 1、用两种不同的顺序计算10000 21n n -=∑，分析其误差的变化 2、已知连分数() 1 01223//(.../)n n a f b b a b a a b =+ +++，利用下面的算法计算f ： 1 1 ,i n n i i i a d b d b d ++==+ (1,2,...,0 i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分 1 041 n n x y dx x =+? (0,1,...,1 n = 4、设2 2 11N N j S j == -∑ ，已知其精确值为1311221N N ?? -- ?+?? （1）编制按从大到小的顺序计算N S 的程序 （2）编制按从小到大的顺序计算N S 的程序 （3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数 三、实验步骤、程序设计、实验结果及分析 1、用两种不同的顺序计算10000 2 1n n -=∑，分析其误差的变化 （1）实验步骤： 分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h （2）程序设计： a.顺序计算 #include 数值计算第一次大作业 实验目的 以Hilbert 矩阵为例，研究处理病态问题可能遇到的困难。 内容 Hilbert 矩阵的定义是 ,() 11/21/31/1/21/31/41/(1)1/31/41/51/(2)1/1/(1)1/(2)1/(21)n i j H h n n n n n n n =????+????=+??????++-?? 它是一个对称正定矩阵，而且()n cond H 随着n 的增加迅速增加，其逆矩阵1,()n i j H α-=， 这里 ,2(1)(1)!(1)!(1)[(1)!(1)!]()!()! i j i j n i n j i j i j n i n j α+-+-+-=+----- 1) 画出ln(())~n cond H n 之间的曲线（可以用任何的一种范数）。你能猜出 ln(())~n cond H n 之间有何种关系吗？提出你的猜想并想法验证。 用行范数 for n=1:50 for i=1:n for j=1:n A(i,j)=1/(i+j-1); B(i,j)=factorial(n+i-1)*factorial(n+j-1)/((i+j-1)*(factorial(i-1)*facto rial(j-1))^2*factorial(n-i)*factorial(n-j)); end end result1=0; for j=1:n result1=result1+A(1,j); end result1=log(result1); result2=0; for i=1:n for j=1:n result2=B(i,j)+result2; end result(i)=log(result2); end m=max(result); 《计算方法》上机指导书 实验1 MATLAB 基本命令 1．掌握MATLAB 的程序设计 实验内容：对以下问题，编写M 文件。 (1) 生成一个5×5矩阵，编程求其最大值及其所处的位置。 (2) 编程求∑=20 1!n n 。 (3) 一球从100米高度自由落下，每次落地后反跳回原高度的一半，再落下。求它在 第10次落地时，共经过多少米？第10次反弹有多高？ 2．掌握MATLAB 的绘图命令 实验内容：对于自变量x 的取值属于[0，3π]，在同一图形窗口画出如下图形。 （1）1sin()cos()y x x =?； （2）21 2sin()cos()3 y x x =-； 实验2 插值方法与数值积分 1. 研究人口数据的插值与预测 实验内容：下表给出了从1940年到1990年的美国人口，用插值方法推测1930年、1965年、2010年人口的近似值。 美国人口数据 1930年美国的人口大约是123,203千人，你认为你得到的1965年和2010年的人口数字精确度如何？ 2．最小二乘法拟合经验公式 实验内容：某类疾病发病率为y ‰和年龄段x (每五年为一段，例如0~5岁为第一段，6~10岁为第二段……)之间有形如bx ae y =的经验关系，观测得到的数据表如下 （1）用最小二乘法确定模型bx ae y =中的参数a 和b 。 （2）利用MATLAB 画出离散数据及拟合函数bx ae y =图形。 3. 复化求积公式 实验内容：对于定积分? +=1 02 4dx x x I 。 （1）分别取利用复化梯形公式计算，并与真值比较。再画出计算误差与n 之间的曲线。 （2）取[0，1]上的9个点，分别用复化梯形公式和复化辛普森公式计算，并比较精度。 题目利用数值计算方法求取基尼系数 姓名与学号 指导教师 年级与专业 所在学院 一、问题综述： 基尼系数（Gini coefficient），是20世纪初意大利学者科拉多·吉尼根据劳伦茨曲线所定义的判断收入分配公平程度的指标。是比例数值，在0和1之间。基尼指数（Gini index）是指基尼系数乘100倍作百分比表示。在民众收入中，如基尼系数最大为“1”，最小等于“0”。前者表示居民之间的收入分配绝对不平均（即所有收入都集中在一个人手里，其余的国民没有收入），而后者则表示居民之间的收入分配绝对平均，即人与人之间收入绝对平等，但这两种情况只出现在理论上；因此，基尼系数的实际数值只能介于0～1之间，基尼系数越小收入分配越平均，基尼系数越大收入分配越不平均。 设右图中的 实际收入分配曲线 （红线）和收入分 配绝对平等线（绿 线）之间的面积为 A，和收入分配绝 对不平等线（蓝 线）之间的面积为 B，则表示收入与 人口之间的比例的基尼系数为 A A+B 。 如果A为零，即基尼系数为0，表示收入分配完全平等（红线和绿线重叠）；如果B为零，则系数为1，收入分配绝对不平等（红线和蓝线重叠）。该系数可在0和1之间取任何值。实际上，一般国家的收入分配，既不是完全平等，也不是完全不平等，而是在两者之间，劳伦茨曲线为一条凸向横轴的曲线。收入分配越趋向平等，劳伦茨曲线的弧度越小（斜度越倾向45度），基尼系数也越小；反之，收入分配越趋向不平等，劳伦茨曲线的弧度越大，那么基尼系数也越大。 基尼系数的调节需要国家通过财政政策进行国民收入的二次分配，例如对民众的财政公共服务支出和税收等，从而让收入均等化，令基尼系数缩小。 基尼系数由于给出了反映居民之间贫富差异程度的数量界线，可以较客观、直观地反映和监测居民之间的贫富差距，预报、预警和防止居民之间出现贫富两极分化。因此得到世界各国的广泛认同和普遍采用。 联合国有关组织规定： ●若低于0.2表示收入平均； ●0.2-0.3表示相对平均； ●0.3-0.4表示相对合理； ●0.4-0.5表示收入差距大； ●0.6以上表示收入差距悬殊。 2013年1月18日，中国国家统计局一次性公布了自2003年以来十年的全国基尼系数。大陆统计局局长马建堂称，按照国际新的统计口径，大陆居民收入的基尼系数，2003年是0.479，2004年是0.473，2005年为0.485，2006年为0.487，2007年为0.484，2008年为0.491，2009年为0.490，2010年为 0.481，2011年为0.477，到2012年的数据是0.474，为2005年以来最低水平，而自2008年起，基尼系数也在逐年下降。而此前西南财大调查数据显示，中国的2012年的基尼系数为0.61，但无论是民间统计的数据还是官方统计的数据，结果都遭到学术界质疑，仍具有争议性。 本文将根据网络上国家统计局的数据，利用上面给出的公式来计算我国从2002年以来的城镇居民基尼系数，并将计算出的数据与现有数据进行比较。 全球基尼系数 计算方法大作业 题目：非线性方程求根的新方法 班级：xxx 学号：xxx 姓名：xxx 非线性方程求根的新方法 一、问题引入 在计算和实际问题中经常遇到如下非线性问题的求解： F(x)=0 (1) 我们经常采用的方法是经典迭代法： 经典迭代方法 不动点迭代方法是一种应用广泛的方法，其加速方法较多，如Stiffensen加速方法的局部收敛阶(以下简称为收敛阶)为2阶；牛顿迭代方法的收敛阶亦为2阶，且与其相联系的一些方法如简化牛顿法、牛顿下山法、弦截法的收敛阶阶数介于1和2之间；而密勒法的收敛阶与牛顿法接近，但计算量较大且涉及零点的选择问题，同时收敛阶也不够理想。 因此本文介绍一种新的迭代方法 从代数角度看，牛顿法和密勒法分别是将f(x)在xk附近近似为一线性函数和二次抛物插值函数，一种很自然的想法就是能否利用Taylor展开，将f(x)在xk附近近似为其他的二次函数?答案是肯定的.其中的一种方法是将f(x)在Xk处展开3项，此时收敛阶应高于牛顿法，这正是本文的出发点． 二、算法推导 设函数f(x)在xk附近具有二阶连续导数，则可将f(x)在xk处进行二阶Taylor展开，方程(1) 可近似为如下二次方程： f(xk)+f’(xk)(x-xk)+2^(-1)f’’(xk)(x-xk)^2=0，(2) 即 2^(-1)f’’(xk)x^2+(f’(xk)-xkf’’(xk))x+2^(-1)f’’(xk)xk^2-xkf’(xk)+f(xk)=0(3) 利用求根公式可得 X=xk-(f’’(xk))^(-1)(f’(xk))-sqrt((f’(xk)^2±2f’’(xk)f(xk)))(4) 其中±符号的选取视具体问题而定，从而可构造迭代公式 X k+1=xk-(f’’(xk))^(-1)(f’(xk))-sqrt((f’(xk)^2±2f’’(xk)f(xk)))(5) 确定了根号前正负号的迭代公式(5)，可称为基于牛顿法和Taylor展开的方法，简记为BNT 方法． 为描述方便起见，以下将f(xk)，f’(xk)，f’’(xk)分别记为f,f’,f’’．首先，二次方程(3)对应于一条抛物曲线，其开口方向由f’’(xk)，x∈U(xk)的符号确定，其中U(xk)为xk的某邻域，其顶点为 P(xk-(f’’)^(-1)f’，fk-(2f’’)^(-1)(f’)^2)．为使(5)式唯一确定x k+1，须讨论根式前正负号的取舍问题．下面从该方法的几何意义分析(5)式中正负号的取舍． 1)当f(xk)=o时，z。即为所求的根． 2)当f(xk)>O时，根据y=f(x)的如下4种不同情形(见图1)确定(5)式中根号前的符号． (a)当f’’(xk) 计算方法大作业利用雅可比迭代求值 班级：021011 学号：02101022 姓名：白文仓 利用雅可比迭代法解方程组 1231231 2352828321361x x x x x x x x x ++=?? +-=??--=? 一、雅可比迭代原理 线性方程组为了计算方便可以用矩阵表示，Ax=b 。一般当A 为低阶稠密矩阵时，用主元消去法解此方程组是有效方法。但是，对于由工程技术中产生的大型稀疏矩阵方程组，利用迭代法求解此方程组比较合适。雅克比迭代法就是众多迭代法中比较早且较简单的一种。 首先将方程组中的系数矩阵A 分解成三部分，即：A = L+D+U ，其中D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。 211 2 00 0n n a L a a ?? ?? ??=???? ??? 11 22 nn a a D a ?? ??? ?=?????? ? 1212000n n a a a U ??????=??????? 之后确定迭代格式，X^(k+1) = B *X^(k) +f ，（这里^表示的是上标，括 号内数字即迭代次数），其中B 称为迭代矩阵，雅克比迭代法中一般记为J 。（k = 0,1，......）。 再选取初始迭代向量X^(0)，开始逐次迭代。 二、收敛性 设A x = b ，其中A=D+L+U 为非奇异矩阵，且对角阵D 也非奇异，则当迭代矩阵J 的ρ（J ）<1时，雅克比迭代法收敛。 三、优缺点 雅克比迭代法的优点明显，计算公式简单，每迭代一次只需计算一次矩 阵和向量的乘法，且计算过程中原始矩阵A 始终不变，比较容易并行计算。然而这种迭代方式收敛速度较慢，而且占据的存储空间较大。 四、算法 设方程组A x = b 的系数矩阵的对角元素()1101,2,,a i n ≠= ，M 为迭代次数容许的最大值， ε 为容许误差。 1），取初始向量 () (0)(0)(0)1 2 ,,,T n x x x x = ，令k=0. 2），对i=1,2，……，n 计算 (1) ()1,1() n k k i i ij j j j i ii x b a x a +=≠=-∑ 3），如果(1)()1 n k k i i i x x ε+=-<∑，则输出(1)k x +，结束；否则执行 4） 4），如果 k M ≥，则不收敛，终止程序；否则 1k k ←+，转2）。 五、 程序 function Jacobi(A,b,max,eps)%max 为迭代次数容许的最大值，eps 为容许误差 A=[5 2 1;2 8 -3;1 -3 -6]; b=[8;21;1]; max=100; eps=0.000001; n=length(A); x=zeros(n,1); x1=zeros(n,1); k=0; while 1 x1(1)=(b(1)-A(1,2:n)*x(2:n,1))/A(1,1); for i=2:n-1 x1(i)=(b(i)-A(i,1:i-1)*x(1:i-1,1)-A(i,i+1:n)*x(i+1:n,1))/A(i,i); end x1(n)=(b(n)-A(n,1:n-1)*x(1:n-1,1))/A(n,n); k=k+1; if sum(abs(x1-x)) 数值分析大作业（2013年5月） 金洋洋（12721512），机自系 1.下列各数都是经过四舍五入得到的近似值，试分别指出它 们的绝对误差限, 相对误差限和有效数字的位数。 X1 =5.420, x 2 =0.5420, x 3=0.00542, x 4 =6000, x 5=50.610? 解：根据定义：如果*x 的绝对误差限 不超过x 的某个数位的半个单位，则从*x 的首位非零数字到该位都是有效数字。 显然根据四舍五入原则得到的近视值，全部都是有效数字。 因而在这里有：n1=4, n2=4, n3=3, n4=4, n5=1 （n 表示x 有效数字的位数） 对x1：有a1=5, m1=1 （其中a1表示x 的首位非零数字，m1表示x1的整数位数） 所以有绝对误差限 143 11 (1)101022 x ε--≤ ?=? 相对误差限 31() 0.510(1)0.00923%5.4201 r x x x εε-?= == 对x2：有a2=5, m2=0 所以有绝对误差限 044 11 (2)101022 x ε--≤ ?=? 相对误差限 42() 0.510(2)0.00923%0.54202 r x x x εε-?= == 对x3：有a3=5, m3=-2 所以有绝对误差限 235 11 (3)101022 x ε---≤ ?=? 相对误差限 53() 0.510(3)0.0923%0.005423 r x x x εε-?= == 对x4：有a4=0, m4=4 所以有绝对误差限 4411(4)1022 x ε-≤?= 相对误差限 4() 0.5 (4)0.0083%6000 4 r x x x εε= = = 对x5：有a5=6, m5=5 所以有绝对误差限 514 11(5)101022 x ε-≤ ?=? 相对误差限 45() 0.510(5)8.3%600005 r x x x εε?= == 工程计算方法与软件应用 本科生大作业 考核方式：考查（成绩按各软件的课外作业成绩综合给出）。 各软件讲完后1~2星期内上交作业。 一、CAD/CAE软件作业（每个学生完成下列任意一题） 题目一： 一端固定支撑，一端集中力的梁，横截面为10x10cm，长为150cm，受集中载荷作用，P=50N。弹性模量E=70GPa，泊松比r=0.2。用ABAQUS 软件建模并计算最大应力和最大位移的位置和大小。 （1）二维；（2）三维 图1梁受力简图 题目二： 图中所示为一个连接件，一端焊接到设备母体上，一端在圆柱销子作用下的圆孔，圆孔下半周受到30 kN的均布载荷作用，用ABAQUS 软件建模并计算最大应力和最大位移的位置和大小。 图2 连接件受力简图 题目三： 如图3所示为一薄壁圆筒，在圆筒中心受集中力F作用，对此进行受力分析，并给出应力、位移云图，并求A、B两点位移。 圆筒几何参数：长度L=0.2m；半径R=0.05m壁厚t=2.5mm。 材料参数：弹性模量E=120Gpa；泊松比0.3 载荷：F=1.5kN。 图3薄壁管受力简图 题目四： 如图4所示为一燃气输送管道截面及受力见图，试分析管道在内部压力作用下的应力场。 几何参数：外径0.6m，内径0.4m，壁厚0.2m 材料参数：弹性模量E=120Gpa；泊松比0.26 载荷P=1Mpa。 图4燃气管受力简图 题目五： 如图5为一三角桁架受力简图，途中各杆件通过铰链链接，杆件材料及几何参数见表1和表2所示，桁架受集中力F1=5kN、F2=2.5kN 作用，求桁架各点位移及反作用力。 图5 三角桁架受力简图 表1 杆件材料参数 表2 杆件几何参数 第四章 问题一 一、问题综述 在离地球表面高度为y处的重力加速度如下： 计算高度y=55000m处的重力加速度值。 二、问题分析 以高度y作为自变量，重力加速度的值为因变量。得到以下信息： f(0)=9.8100; f(30000)=9.7487; f(60000)=9.6879; f(90000)=9.6278; f(120000)=9.5682; 本题要求的就是f(55000)的值。 以下将采用课堂中学到的Lagrange插值多项式法、Newton插值多项式法、分段低次插值法和样条插值法求解该问题。 三、问题解决 1. lagrange插值多项式法 对某个多项式函数，已知有给定的k+ 1个取值点： 其中对应着自变量的位置，而对应着函数在这个位置的取值。 假设任意两个不同的x j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为： 其中每个为拉格朗日基本多项式（或称插值基函数），其表达式为： 拉格朗日基本多项式的特点是在上取值为1，在其它的点上取值为0。 源程序lagrange.m function [c,f]=lagrange(x,y,a) % 输入：x是自变量的矩阵；y是因变量的矩阵；a是要计算的值的自变量； % 输出：c是插值多项式系数矩阵；f是所求自变量对应的因变量； m=length(x); l=zeros(m,m); % l是权矩阵 f=0; for i=1:m v=1; for j=1:m if i~=j v=conv(v,poly(x(j)))/(x(i)-x(j)); % v是l_i(x)的系数矩阵 end end l(i,:)=v; % l矩阵的每一行都是x从高次到低次的系数矩阵 end c=vpa(y*l,10); % 对应阶次的系数相加，乘以y，显示10位有效数字 for k=1:m f=f+c(k)*a^(m-k); end 输入矩阵 x=[0 30000 60000 90000 120000] y=[9.81 9.7487 9.6879 9.6278 9.5682] a=55000 再运行源函数，可得： c = [ -2.057613169e-23, 4.938271605e-18, -3.703703702e-14, -0.000002046111111, 9.81] f = 9.6979851723251649906109417384537数值计算方法大作业
数值分析大作业三 四 五 六 七
计算方法上机作业
西工大计算方法作业答案
西安交通大学计算方法B大作业资料
2020年奥鹏吉大网络教育《计算方法》大作业解答
数值分析大作业
西安交大计算方法b大作业课件
计算方法上机实习题大作业(实验报告).
计算方法大作业第一次
计算方法作业2
数值计算方法大作业
计算方法大作业非线性方程求根的新方法
计算方法大作业----------雅可比
上海大学_王培康_数值分析大作业
工程计算方法及软件应用--本科生考查大作业
数值计算方法第4次作业
- 2020年奥鹏吉大网络教育《计算方法》大作业解答
- 2018-2019学年第二学期期末考试《计算方法》大作业参考答案
- 数值计算方法大作业
- 西安交大计算方法b大作业
- 2019-2020学年第一学期期末考试《计算方法》大作业答案
- 上海交通大学计算方法大作业
- 2019-2020学年第一学期期末考试《计算方法》大作业答案
- 计算方法大作业1 克服Runge现象
- 北京科技大学计算方法大作业
- 计算方法大作业
- 计算方法
- 西工大计算方法作业答案
- 北京理工大学数值计算方法大作业数值实验1
- 西交计算方法A上机大作业
- 浙江大学计算方法大作业第九题
- 计算方法大作业
- 计算方法各章作业答案
- 计算方法上机实习题大作业(实验报告).
- 浙江大学计算方法大作业第八题
- 计算方法大作业作业((北京科技大学研究生结课考试)