第一篇:数值分析教教案18
4.1.3 Newton法
1.算法的基本思想
图4-5 牛顿法原理示意图
将非线性方程线性化,以线性方程的解逐步逼近非线性方程的解,这就是Newton法的基本思想。
把函数f(x)在某一初始值x0点附近展开成Taylor级数有
f(x0)f(x)f(x0)(xx0)f(x0)(xx0)取其
2!线性部分,近似地代替函数f(x)可得方程的近似式:
2f(x)f(x0)(xx0)f(x0)0
设f(x0)0,解该近似方程可得:
f(x0)x1x0f(x0)
可作为方程(4-1)的近似解。重复以上过程,得迭代公式
xk1f(xk)xk(k0,1,2,)(4-8)
f(xk)按式(4-8)求方程(4-1)的近似解称为Newton法。Newton法也是一种不动点迭代,其迭代函数为
f(x)g(x)xf(x)
如图(4-5)所示,从几何上看,yf(x0)f(x0)(xx0)为曲线yf(x)过点(x0,f(x0))的切线,x1为切线与x轴交点,x2则
x轴的交点。如此继续下去,xk1为曲线上点(x,f(x))处的切线与x轴的交点。因此Newton法是以曲线的切线与x轴的交点作为曲线与x轴的交点的近似,故是曲线上点(x1,f(x1))处的切线与
kkNewton法又称为切线法。
2.切线法的收敛性
理论可以证明,在有根区间[a,b]上,f(x)0,f(x)连续且不变号,则只要选取的初始近似根
x0满足f(x0)f(x0)0,切线法必定收敛。它的收敛速度可如下推出。
方程f(x)0可以等价地写成f(x)(xx)f(x),若f(x)0f(x)f(x)移项可得xxf(x)。设g(x)xf(x),两边求导得f(x)f(x)g(x)f(x)0g(x)0。xx,则必得 2,代入[f(x)]另一方面,比较迭代公式xk1f(xk)xkf(xk)和f(x)g(x)xg(x)xg(x)x可知。把函数在点展k1kf(x)成泰勒级数,只取到二阶导数则有:
g(x)2g(xk)g(x)g(x)(xkx)(xkx)由
2g(x)g(xk)g(x)(xkx)2,移
2xk1g(x)0,所以有xk1于xg(x),得出 项并注意
g(x)xk1g(x)xk1x(xkx)2
2为了将式中的g(x)换成ff(x)f(x)(x),对g(x)[f(x)]2两边求导,并代入f(x)0,则有:
f(x)g(x)f(x)*将它代入前式得出
f(x)2xk1x(xkx)(4-9)2f(x)f(x)2f(x)是个常数,式(4-9)表明用牛顿迭代公式在某次算得的误差,与上次误差的平方成正比,可见牛顿迭代公式的收敛速度很快,但计算实践表明,当初值不够好时,Newton法可能发散。一般可由问题的实际背景来预测或由对分区间法求得较好的初始值。
3.Newton迭代公式Matlab实现
按照算法4-4编写迭代法的Matlab程序(函数名:Newton.m).function [p1,err,k,y]=newton(f,df,p0,delta,max1)% f是非线性函数 % df是f的微商 % p0是初始值
% delta是给定允许误差 % max1是迭代的最大次数
% p1是牛顿法求得的方程的近似值 % err是p0的误差估计 % k是迭代次数 % y=f(p1)p0,feval('f',p0)for k=1:max1 p1=p0-feval('f',p0)/feval('df',p0);err=abs(p1-p0);p0=p1;p1,err,k,y=feval('f',p1)if(err 3x3x2的近似值,给定一个初始值【例4-3】求方程p01.2,误差界为106。 首先,在MATLAB命令窗口输入: fplot('[x^3-3*x+2,0]',[-2.5 2.5]);grid;回车得到如图4-6所示图形,即可知函数f(x)与x轴有交点,也就是说有根,并且从图中能够大致估算到根的位置。 3f(x)x3x2: 先用一个名为f.m的文件定义函数function y=f(x)y=x^3-3*x+2; 2df(x)3x3: 再用一个名为df.m的文件定义函数的微商function y=df(x)y=3*x^2-3;然后在MATLAB命令窗口输入: >> newton('f','df',1.2,10^(-6),10) 图4-6 f(x)与x轴交点显示图 回车得到如下结果(只给出了最后两次迭代结果): …… k = 9 p1 = 1.0004 err = 4.1596e-004 k = 10 p1 = 1.0002 err = 2.0802e-004 ans = 1.0002 4.1.4弦截法 1.弦截法基本原理 图4-7弦截法原理示意图 f(xk)的值,当f比较复杂时难以实现,下面要介绍的弦截法用f(x)在两个点上的值构造一次插值函用牛顿法解非性方程要知道数来回避微商的计算。 给定非线性方程f(x)0,选定曲线 yf(x)上的两个点P0(x0,f(x0)),P1(x1,f(x1)),过这两个点作一条直线P0P1,则直f(x1)f(x0)(xx1)。当线方程yf(x1)x1x0f(x0)f(x1)时,f(x1)(x1x0)x2x1xx直线P0P与轴交点为,这时用2作为1f(x1)f(x0)曲线yf(x)与x轴交点的近似值,显然: x2xmin(x1x,x0x) 这里***x*为f(x)0的精确解。然后用 P1(x1,f(x1)),P2(x2,f(x2)),构造直线P1P2,重复上述步骤,可以求出x3。 如此进行下去,得到迭代格式: f(xk)xk1xk(xkxk1)(k0,1,)(4-10)f(xk)f(xk1)f(xk)f(xk1)迭代格式(4-10)实际上就是用差商取代牛顿公xkxk1式xk1f(xk)xk(x)的结果,所以弦截法可以看成f中的微商f(xk)牛顿法的一种变形。 弦截法与牛顿虽然都是非线性方程线性化的方法,但牛顿法在计算xk1时只用到xk的值,但是弦截法要用到xk和xk1的值,故这就要给定两个初值。切线法具有恒收敛且收敛速度快的优点,但需要求出函数的导数;弦截法不需要求导数,收敛速度很快,但是需要知道两个近似的初始值才能作出弦,要求的初始值条件较多。2.弦截法的MATLAB实现 function [p1,err,k,y]=secant(f,p0,p1,delta,max1)% f是给定的非线性函数 % p0,p1为初始值 % delta为给定误差界 % max1是迭代次数的上限 % p1为所求得的方程的近似解 % err为p1-p0绝对值 % k为所需要的迭代次数 % y=f(p1)k=0,p0,p1,feval('f',p0),feval('f',p1)for k=1:max1 p2=p1-feval('f',p1)*(p1-p0)/(feval('f',p1)-feval('f',p0));err=abs(p2-p1);p0=p1;p1=p2;k,p1,err, y=feval('f',p1)if(err xx20,给定初值为 3p01.5,p11.52,误差界为106。 首先,在MATLAB命令窗口输入: fplot('[x^3-x+2,0]',[-2.5 2.5]);grid;回车得到如图4-8所示图形,即可知函数f(x)与x轴有交点,也就是说有根,并且从图中能够大致估算到根的位置。 3f(x)xx2; 解:先用一个名为f.m的文件定义function y=f(x)y=x^3-x+2;然后在MATLAB命令窗口输入: >> secant('f',-1.5,-1.52,10^(-6),11)k =0 p0 =-1.5000 p1 =-1.5200 ans =0.1250 ans =0.0082 k =1 p1 =-1.5214 err =0.0014 y =-1.3633e-004 k =2 p1 =-1.5214 err =2.2961e-005 y =1.4454e-007 k =3 p1 =-1.5214 err =2.4318e-008 y =2.5460e-012 ans =-1.5214 这就表明经过了3次迭代得到满足精足要求的近似解xx31.5214,且f(x3)2.5461012。 图4-8 f(x)与x轴交点显示图 4.1.5抛物线法 1.抛物法的基本原理 设已知方程f(x)=0的三个近似根xk,xk-1,xk-2,我们以这三点为节点构造二次插值多项式P(x)的一个零点xk+1作为新的近似根,这样确定的迭代过程称抛物线法。 在几何图形上,这种方法的基本思想是用抛物线与 x轴的交点xk+1作为所求根的近似位置(图6-7)。 图4-9抛物法图形 现在推导抛物线法的计算公式.插值多项式 P(x)=f(xk)+ f[xk,x k-1](x-xk)+f[xk,xk-1,xk-2](x-xk)(x-xk-1) 有两个零点: xk+1=xk-式中 2f(xk)ω±ω2-4f(xk)f[xk,xk-1,xk-2](4-11) ω=f[xk,xk-1]+f[xk,xk-1,xk-2](xk-xk-1) 为了从(4-11)式定出一个值取舍问题. 在xk+1,我们需要讨论根式前正负号的xk,xk-1,xk-2三个近似根中,自然假定xk更接近所求的根x*,这时,为了保证精度,我们选定(4-11)中较接近xk的一个值作为新的近似根xk+1,为此,只要取根式前的符号与ω的符号相同。 2.抛物法的计算步骤 给定非线性方程f(x)0,误差界ε,迭代次数上限N。(1)计算ω=(2)计算f[xk,xk- 1]+f[xk,xk- 1,xk- 2](xk-xk- 1) 2f(xk),代入: ω±ω2-4f(xk)f[xk,xk-1,xk-2]xk+ 1=xk-2f(xk)ω±ω2-4f(xk)f[xk,xk-1,xk-2] 得出xk+1的值后,计算f(xk+1)。 *x-xx(3)若k+1≈xk+1;否则的话,令: k≤ε则迭代停止,取(xk- 2,xk- 1,xk,f(xk- 2),f(xk- 1),f(xk))=(xk- 1,xk,xk+ 1,f(xk- 1),f(xk),f(xk+ 1)),n=n+1。 (4)如果迭代次数k>N,则认为该迭代格式对于所选初值不收敛,迭代停止;否则的话重返步骤(2)。3.抛物法Matlab实现 按照算法编写迭代法的Matlab程序 function root=Parabola(f,a,b,x,eps)% f是非线性函数 % a为有根区间的左限 % b为有根区间的右限 % eps为根的精度 % root为求出的函数零点 % x为初始迭代点的值 if(nargin==4)eps=1.0e-4;end f1=subs(sym(f),findsym(sym(f)),a);% subs命令是符号函数赋值 f2=subs(sym(f),findsym(sym(f)),b);% findsym确定自由变量是对整个矩阵进行的 if(f1==0)root=a;end if(f2==0)root=b;end if(f1*f2>0)disp('两端点函数值乘积大于0!');return;else tol=1;fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);fx=subs(sym(f),findsym(sym(f)),x);d1=(fb-fa)/(b-a);d2=(fx-fb)/(x-b);d3=(d2-d1)/(x-a);B=d2+d3*(x-b);root=x-2*fx/(B+sign(B)*sqrt(B^2-4*fx*d3));t=zeros(3);t(1)=a;t(2)=b;t(3)=x;while(tol>eps)t(1)=t(2);t(2)=t(3);t(3)=root;f1=subs(sym(f),findsym(sym(f)),t(1));f2=subs(sym(f),findsym(sym(f)),t(2));f3=subs(sym(f),findsym(sym(f)),t(3));d1=(f2-f1)/(t(2)-t(1));d2=(f3-f2)/(t(3)-t(2));d3=(d2-d1)/(t(3)-t(1));B=d2+d3*(t(3)-t(2));root=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3));tol=abs(root-t(3));end end 【例4-5】用抛物法求解方程lgx+x=2在区间[1,4]内的一个根。 首先,在MATLAB命令窗口输入: fplot('[sqrt(x)+log(x)-2]',[1 4]);grid;回车得到如图4-10所示图形,即可知函数f(x)与x轴有交点,也就是说有根,并且从图中能够大致估算到根的位置。 1.510.50-0.5-1 图4-10 f(x)与x轴交点显示图 然后在MATLAB命令窗口输入: >> r=Parabola(' sqrt(x)+log(x)-2',1,4,2)回车得到如下结果: r = 1.8773 从结果中可以知道,方程lgx+11.522.533.54x=2的一个根x=1.8773。 数值分析学习感想 一个学期的数值分析,在老师的带领下,让我对这门课程有了深刻的理解和感悟。这门课程是一个十分重视算法和原理的学科,同时它能够将人的思维引入数学思考的模式,在处理问题的时候,可以合理适当的提出方案和假设。他的内容贴近实际,像数值分析,数值微分,求解线性方程组的解等,使数学理论更加有实际意义。 数值分析在给我们的知识上,有很大一部分都对我有很大的帮助,让我的生活和学习有了更加方便以及科学的方法。像第一章就讲的误差,在现实生活中,也许没有太过于注意误差,所以对误差的看法有些轻视,但在学习了这一章之后,在老师的讲解下,了解到这些误差看似小,实则影响很大,更如后面所讲的余项,那些差别总是让人很容易就出错,也许在别的地方没有什么,但是在数学领域,一个小的误差,就很容易有不好的后果,而学习了数值分析的内容,很容易就可以将误差锁定在一个很小的范围内,在这一范围内再逼近,得出的近似值要准确的多,而在最开始的计算中,误差越小,对后面的影响越小,这无疑是好的。 数值分析不只在知识上传授了我很多,在思想上也对我有很大的影响,他给了我很多数学思想,很多思考的角度,在看待问题的方面上,多方位的去思考,并从别的例子上举一反 三。像其中所讲的插值法,在先学习了拉格朗日插值法后,对其理解透彻,了解了其中的原理和思想,再学习之后的牛顿插值以及三次样条插值等等,都很容易的融会贯通,很容易的就理解了其中所想,他们的中心思想并没有多大的变化,但是使用的方式却是不同的,这不仅可以学习到其中心内容,还可以去学习他们的思考方式,每个不同的思考方式带来的都是不同的算法。而在看待问题上,不同的思考方式总是可以快速的全方位的去看透彻问题,从而知道如何去解决。 在不断的学习中,知识在不断的获取,能力在不断的提升,同时在老师的不懈讲解下,我逐渐的发现数值分析所涵盖的知识面特别的广泛,而我所需要学习的地方也更加的多,自己的不足也在不断的体现,我知道这只是我刚刚接触到了数学的那一角,在以后我还会接触到更多,而这求知的欲望也在不停的驱赶我,学习的越多,对今后的生活才会有更大的帮助。 计算132 2013014923 张霖篇二:数值分析学习报告 数值分析学习心得报告 班级:11级软工一班 姓名: * * * 学号: 20117610*** 指导老师:* * * 学习数值分析的心得体会 无意中的一次选择,让我接触了数值分析。 作为这学期的选修课,我从内心深处来讲,数值分析真的有点难。感觉它是在高等数学和线性代数的基础上,又加深了探讨。虽然这节课很难,我学的不是很好,但我依然对它比较感兴趣。下面就具体说说我的学习体会,让那些感兴趣的同学有个参考。学习数值分析,我们首先得知道一个软件——matlab。matrix laboratory,即矩阵实验室,是math work公司推出的一套高效率的数值计算和可视化软件。它是当今科学界最具影响力、也是最具活力的软件,它起源于矩阵运算,并高速发展成计算机语言。它的优点是强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面、便捷的与其他程序和语言接口。 根据上网搜集到的资料,你就会发现matlab有许多优点: 首先,编程简单使用方便。到目前为止,我已经学过c语言,机器语言,java语言,这三个语言相比,我感觉c语言还是很简单的一种编程语言。只要入门就很好掌握,但是想学精一门语言可不是那么容易的。惭愧的说,到目前为止,我依然处于入门阶段,只会编写小的简单的程序,但是班里依然还是有学习好的。c语言是简单且容易掌握的,但是,matlab的矩阵和向量操作功能是其他语言无法比拟的。在matlab环境下,数组的操作与数的操作一样简单,基本数据单元是不需要指定维数的,不需要说明数据类型的矩阵,而其数学表达式和运算规则与通常的习惯相同。 其次,函数库可任意扩充。众所周知,c语音有着丰富的函数库,我们可以随时调用,大大方便了程序员的操作。可是作为it人士的你知道吗,由于matlab语言库函数与用户文件的形式相同,用户文件可以像库函数一样随意调用,所以用户可任意扩充库函数。这是不是很方便呢? 接着,语言简单内涵丰富。数值分析所用的语言中,最重要的成分是函数,其一般形式为:function[a,b,c??]=fun(d,e,f??),你也发现了吧,这样的语音是不是很容易掌握呢!fun是自定义的函数名,只要不与库函数想重,并且符合字符串书写规则即可。 然后是丰富的工具箱。由于matlab 的开放性,许多领域的专家都为matlab 编写了各种程序工具箱。这些工具箱提供了用户在特别应用领域所需的许多函数,这使得用户不必花大量的时间编写程序就可以直接调用这些函数,达到事半功倍的效果。不过你得提前知道这些工具箱,并且会使用。 最后,我们来说一下matlab的运算。利用matlab可以做向量与矩阵的运算,与普通加减运算几乎相似。 矩阵乘法用 “ * ” 符号表示,当a矩阵列数与b矩阵的行数相等时,二者可以进行乘法运算,否则是错误的。如果a或b是标量,则a*b返回标量a(或b)乘上矩阵b(或a)的每一个元素所得的矩阵。 对n×m阶矩阵a和p×q阶矩阵b,a和b的kronecher乘法运算可定义为: kronecker乘法的matlab命令为c=kron(a,b):例如,在matlab中输入: a=[1 2;3 4];b=[1 3 2;2 4 6];c=kron(a,b)则程序会给出相应的答案 c = 1 3 2 2 6 4 2 4 6 4 8 12 3 9 6 4 12 8 6 12 18 8 16 24 这就充分的考验了我们的实际动手能力,当然运用一般的计算方法能算出结果,但相对来说没有用它来运算节省时间,其他算法又很不方便。上面介绍了matlab的特点与使用方法,接着我们要说它的程序设计,其实跟c语言相比,它们的程序设计都差不多。 大家都知道,matlab与其它计算机语言一样,也有控制流语句。而控制流语句本身,可使原本简单地在命令行中运行的一系列命令或函数,组合成为一个整体—程序,从而提高效率。以下是具体的几个例子,看过之后,你会发现,matlab的控制流语句跟其他计算机真的很相似: (1)for 循环for循环的通用形式为:for v=expressionstatementsend其中expression 表达式是一个矩阵,因为matlab中都是矩阵,矩阵的列被一个接一个的赋值到变量v,然后statements语句运行。 (2)while 循环while循环的通用形式为:while v=expressionstatementsend当expression的所有运算为非零值时,statements 语句组将被执行。如果判断条件是向量或矩阵的话,可能需要all 或any函数作为判断条件。(3)if和break语句通用形式为:if 条件1,命令组1;elesif条件2,命令组2;??;else命令组k;endbreak%中断执行,用在循环语句内表示跳出循环。对于数值分析这节课,我的理解是:只要学习并掌握好matlab,你就已经成功了。因此说,matlab是数学分析的基础。另外,自我感觉这是一个很好的软件,其语言简便,实用性强。但是作为一个做新手,想要学习好这门语言,还是比较困难的。在平常的上机课中,虽然我没有问过老师,但是我向那些学习不错的学生还是交流了许多,比如说,张**,贾**,还有那个皮肤白白的女生。跟他们交流,我确实学到不少有用的东西。但是,毕竟没有他们学得好,总之,在我接触这门语言的这些天,除了会画几个简单的三维图形,其他的还是有待提高。在这个软件中,虽然有help,但大家不要以为有了这个就万事大吉了,反而,从另一个方面也对我们大学生提出了两个要求——充实的课外基础和良好的英语基础。在现代,几乎所有好的软件都是来自国外,假如你不会外语,想学好是非常难的,即使高考中的英语比重降低了,但我们依旧得学好。这样我们才能走得更远。 其实想要学习好一们语言,不能只靠老师,靠朋友,关键是自己。每个人内心深处都是有抵触意识的,不可能把老师的所有都学到。其实,我发现学习数值分析这门课,不光是学习一种语言,一些知识,更重要的是学习一种方法,一种学习软件的方法,还有学习的态度。 在最后,我想说的是,谢谢郭老师的辛勤付出,我们每个学生都会看在眼里记在心里的,谢谢您。篇三:数值分析学习总结感想 数值分析学习感想 摘要:数值分析主要介绍现代科学计算中常用的数值计算方法及其基本原理,研究并解决数值问题的近似解,是数学理论与计算机和实际问题的有机结合。随着科学技术迅速发展,运用数学方法解决工程技术领域中的实际问题,已经得到普遍重视。 作为这学期的考试课,在我最初接触这门课时,我感到了很困难,因为无论是高数还是线性代数我都放下了很久,而我感觉数值分析是在高等数学和线性代数的基础上,又加深了探讨。虽然这节课很难,但是在老师不断地引导和讲授下,我逐渐对其产生了兴趣。在老师的反复讲解下,我发现我被它吸引了,因为它不仅是单纯的学科,还教会了我许多做人生活的道理。 首先,数值分析这门课程是一个十分重视算法和原理的学科,同时它能够将人的思维引入数学思考的模式,在处理问题的时候,可以合理适当的提出方案和假设。他的内容贴近实际,像数值分析,数值微分,求解线性方程组的解等,使数学理论更加有实际意义。 数值分析在给我们的知识上,有很大一部分都对我有很大的帮助,让我的生活和学习有了更加方便以及科学的方法。像第一章就讲的误差,在现实生活中,也许没有太过于注意误差,所以对误差的看法有些轻视,但在学习了这一章之后,在老师的讲解下,了解到这些误差看似小,实则影响很大,更如后面所讲的余项,那些差别总是让人很容易就出错,也许在别的地方没有什么,但是在数学领域,一个小的误差,就会有很大的差别,而学习了数值分析的内容,很容易就可以将误差锁定在一个很小的范围内,在这一范围内再逼近,得出的近似值要准确的多,而在最开始的计算中,误差越小,对后面的影响越小,这无疑是好的。数值分析中,“以点带面”的思想也深深影响了我。这里的“点”是根本,是主线。在第二章学习插值法的时候是以拉格朗日插值、牛顿插值为主线,然后逐渐展开介绍艾尔米特插值、分段低次插值和三次样条插值。在学习中只要将研究拉格朗日插值和牛顿插值的基本原理、基本方法理解透彻,其他的插值方法就基本掌握了。第四章处理数值积分和数值微分的基本方法是逼近法,只要将函数逼近的基本思想理解好,掌握起来就会得心应手;第六第七章是以迭代法为主线来求解线性方程组和非线性方程组的。在学习过程组只要将迭代法的相关原理掌 握好,便能掌握第六第七章。总的来数,数值分析所涉及到数学中很多学科的知识,内容比较复杂,因此在学习过程中一定要将基本原理、基本算法理解透,然后再逐步推广。同样在生活中每件事情都有它的主线,只要抓住这条主线再难的事情也会迎刃而解。 还比如“等价转化”的思想,这里的“等价”不是完全意义上的“等价”,是指在转化前后转化的主体主要特征值没有变。插值法的思想就是抓住已知函数或者已知点的几个主要特征,用另一个具备主要特征的简单函数来代替原函数或拟合已知数据点。实际生活中也有很多类似情况,已知事件或者面临的情况往往是复杂的,常常不能直接用数学方法直接研究,我们可以做的就是抓住已经事件的主要特征转化为数学模型来建立。 在不断的学习中,知识在不断的获取,能力在不断的提升,同时在老师的耐心讲解下,我逐渐的发现数值分析所涵盖的知识面特别的广泛,而我所需要学习的地方也更加的多,自己的不足也在不断的体现,我知道这只是我刚刚接触到了数学的那一角,在以后我还会接触到更多,而这求知的欲望也在不停的驱赶我,学习的越多,对今后的生活才会有更大的帮助。 希望在将来,通过反复的实践能加深我的理解,在明年的这个时候我能有更多的感悟。同时,因为十五周的学习时间太短加上我的基础薄弱,我决定明年继续来旁听老师的课程,达到进一步学习,加深理解的目的。 数值分析课程论文: 数值分析学习心得感悟 姓名:崔俊毅 学号:2015210211 专业:防灾减灾专硕 院系:土木工程学院篇四:数值分析学习报告 数值分析学习心得报告 班级:姓名: 学号: ************ *** *********** 学习数值分析的心得体会 数值分析是一门利用计算机求解数学问题数值解的课程,有很强的理论性和实践性,无意中的一次选择,让我接触了数值分析。随着科学技术的发展,提出了大量复杂的数值计算问题,在建立电子计算机成为数值计算的主要工具以后,它以数字计算机求解数学问题的理论和方法为研究对象。有可靠的理论分析,要有数值实验,并对计算的结果进行误差分析。数值分析的主要内容包括插值法,函数逼近,曲线拟和,数值积分,数值微分,解线性方程组的直接方法,解线性方程组的迭代法,非线性方程求根,常微分方程的数值解法。 作为这学期的选修课,我从内心深处来讲,数值分析真的有点难。感觉它是在高等数学和线性代数的基础上,又加深了探讨。虽然这节课很难,我学的不是很好,但我依然对它比较感兴趣。下面就具体说说我的学习体会,让那些感兴趣的同学有个参考。学习数值分析,我们首先得知道一个软件——matlab。matrix laboratory,即矩阵实验室,是math work公司推出的一套高效率的数值计算和可视化软件。它是当今科学界最具影响力、也是最具活力的软件,它起源于矩阵运算,并高速发展成计算机语言。它的优点是强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面、便捷的与其他程序和语言接口。 根据上网搜集到的资料,你就会发现matlab有许多优点: 首先,编程简单使用方便。到目前为止,我已经学过c语言,机器语言,java语言,这三个语言相比,我感觉c语言还是很简单的一种编程语言。只要入门就很好掌握,但是想学精一门语言可不是那么容易的。惭愧的说,到目前为止,我依然处于入门阶段,只会编写小的简单的程序,但是班里依然还是有学习好的。c语言是简单且容易掌握的,但是,matlab的矩阵和向量操作功能是其他语言无法比拟的。在matlab环境下,数组的操作与数的操作一样简单,基本数据单元是不需要指定维数的,不需要说明数据类型的矩阵,而其数学表达式和运 算规则与通常的习惯相同。 其次,函数库可任意扩充。众所周知,c语音有着丰富的函数库,我们可以随时调用,大大方便了程序员的操作。可是作为it人士的你知道吗,由于matlab语言库函数与用户文件的形式相同,用户文件可以像库函数一样随意调用,所以用户可任意扩充库函数。这是不是很方便呢? 接着,语言简单内涵丰富。数值分析所用的语言中,最重要的成分是函数,其一般形式为:function[a,b,c??]=fun(d,e,f??),你也发现了吧,这样的语音是不是很容易掌握呢!fun是自定义的函数名,只要不与库函数想重,并且符合字符串书写规则即可。 然后是丰富的工具箱。由于matlab 的开放性,许多领域的专家都为matlab 编写了各种程序工具箱。这些工具箱提供了用户在特别应用领域所需的许多函数,这使得用户不必花大量的时间编写程序就可以直接调用这些函数,达到事半功倍的效果。不过你得提前知道这些工具箱,并且会使用。 最后,我们来说一下matlab的运算。利用matlab可以做向量与矩阵的运算,与普通加减运算几乎相似。 矩阵乘法用 “ * ” 符号表示,当a矩阵列数与b矩阵的行数相等时,二者可以进行乘法运算,否则是错误的。如果a或b是标量,则a*b返回标量a(或b)乘上矩阵b(或a)的每一个元素所得的矩阵。 对n×m阶矩阵a和p×q阶矩阵b,a和b的kronecher乘法运算可定义为: kronecker乘法的matlab命令为c=kron(a,b):例如,在matlab中输入: a=[1 2;3 4];b=[1 3 2;2 4 6];c=kron(a,b)则程序会给出相应的答案 c = 1 3 2 2 6 4 2 4 6 4 8 12 3 9 6 4 12 8 6 12 18 8 16 24 这就充分的考验了我们的实际动手能力,当然运用一般的计算方法能算出结果,但相对来说没有用它来运算节省时间,其他算法又很不方便。上面介绍了matlab的特点与使用方法,接着我们要说它的程序设计,其实跟c语言相比,它们的程序设计都差不多。 大家都知道,matlab与其它计算机语言一样,也有控制流语句。而控制流语句本身,可使原本简单地在命令行中运行的一系列命令或函数,组合成为一个整体—程序,从而提高效率。以下是具体的几个例子,看过之后,你会发现,matlab的控制流语句跟其他计算机真的很相似: (1)for 循环for循环的通用形式为:for v=expressionstatementsend其中expression 表达式是一个矩阵,因为matlab中都是矩阵,矩阵的列被一个接一个的赋值到变量v,然后statements语句运行。 (2)while 循环while循环的通用形式为:while v=expressionstatementsend当expression的所有运算为非零值时,statements 语句组将被执行。如果判断条件是向量或矩阵的话,可能需要all 或any函数作为判断条件。 (3)if和break语句通用形式为:if 条件1,命令组1;elesif条件2,命令组2;??;else命令组k;endbreak%中断执行,用在循环语句内表示跳出循环。对于数值分析这节课,我的理解是:只要学习并掌握好matlab,你就已经成功了。因此说,matlab是数学分析的基础。另外,自我感觉这是一个很好的软件,其语言简便,实用性强。但是作为一个做新手,想要学习好这门语言,还是比较困难的。在平常的上机课中,虽然我没有问过老师,但是我向那些学习不错的学生还是交流了许多,跟他们交流,我确实学到不少有用的东西。但是,毕竟没有他们学得好,总之,在我接触这门语言的这些天,除了会画几个简单的三维图形,其他的还是有待提高。在这个软件中,虽然有help,但大家不要以为有了这个就万事大吉了,反而,从另一个方面也对我们大学生提出了两个要求——充实的课外基础和良好的英语基础。在现代,几乎所有好的软件都是来自国外,假如你不会外语,想学好是非常难的,即使高考中的英语比重降低了,但我们依旧得学好。这样我们才能走得更远。其实想要学习好一们语言,不能只靠老师,靠朋友,关键是自己。每个人内心深处都是有抵触意识的,不可能把老师的所有都学到。其实,我发现学习数值分析这门课,不光是学习一种语言,一些知识,更重要的是学习一种方法,一种学习软件的方法,还有学习的态度。 数值分析是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科,是数学的一个分支,它以数字计算机求解数学问题的理论和方法为研究对象。在科学研究和工程技术中有许多问题可归结为求解方程组的问题。本文主要讨论了插值法求函数,解线性方程组的求解方法,非线性方程组的解法及微分方程的解法,并通过在电流回路和单晶硅提拉过程中分析应用。进一步体现了数值分析的广泛应用,实际上由于误差的存在,一些问题只能求得近似解。对于良态方程组,只要求解方法稳定,即可得到比较满意的计算结果。但对于病态方程组,即使使用稳定性好的算法求解也未必理想,还需进一步的研究。总之,数值分析可以通过计算方法进行一种比较完善的构造,使之更普遍化,能够有举一反三的思想,能够解决一些实际中难解的问题,应用到各个领域。 在最后,我想说的是,谢谢老师的辛勤付出,我们每个学生都会看在眼里记在心里的,谢谢您。篇五:数值分析期末总结论文,程序界面 数值计算方法论文 论文名称:数值计算方法期末总结 学 号: 姓 名:完成时间: 摘要:数值计算方法是数学的一个重要分支,以用计算机求解数学问题的理论和方法为研究对象。本文是我对本学期数值分析这门课程中所学到的内容以及所作的工作的总结。通过一学期的学习,我深入学习了线性方程组的解法,非线性方程的求根方法,矩阵特征值与特征向量的计算,函数的插值方法,最佳平方逼近,数值积分与数值微分,常微分方程初值问题的数值解法。通过陶老师课堂上的讲解和课下的上机训练,对以上各个章节的算法有了更深刻的体会。最后做了程序的演示界面,使得程序看起来清晰明了,便于查看与修改。通过本学期的学习。 关键词:数值计算方法、演示界面 第一章 前言 随着电子计算机的普及与发展,科学计算已成为现代科学的重要组成部分,因而数值计算方法的内容也愈来愈广泛和丰富。通过本学期的学习,主要掌握了一些数值方法的基本原理、具体算法,并通过编程在计算机上来实现这些算法。 第二章 基本概念 2.1算法 算法是指由基本算术运算及运算顺序的规定构成的完整的解题步骤。算法可以使用框图、算法语言、数学语言、自然语言来进行描述。具有的特征:正确性、有穷性、适用范围广、运算工作量少、使用资源少、逻辑结构简单、便于实现、计算结果可靠。2.2 误差 计算机的计算结果通常是近似的,因此算法必有误差,并且应能估计误差。误差是指近似值与真正值之差。绝对误差是指近似值与真正值之差或差的绝对值;相对误差:是指近似值与真正值之比或比的绝对值。误差来源见表2.1 表 第三章 泛函分析 2.1泛函分析概要 泛函分析(functional analysis)是研究“函数的函数”、函数空间和它们之间变换(映射)的一门较新的数学分支,隶属分析数学。它以各种学科为具体背景,在集合的基础上,把客观世界中的研究对象抽象为元素和空间。如:距离空间,赋范线性空间,内积空间。2.2 范数 范数,是具有“长度”概念的函数。在线性代数、泛函分析及相关的数学领 域,泛函是一个函数,其为矢量空间内的所有矢量赋予非零的正长度或大小。 这里以cn空间为例,rn空间类似。最常用的范数就是p-范数。若,那么 当p取1,2,∞的时候分别是以下几种最简单的情形: 1-范数:║x║1=│x1│+│x2│+?+│xn│ 2-范数:║x║2=(│x1│2+│x2│2+?+│xn│2)1/2 ∞-范数:║x║∞=max(│x1│,│x2│,?,│xn│) 其中2-范数就是通常意义下的距离。 对于这些范数有以下不等式:║x║∞ ≤ ║x║2 ≤ ║x║1 ≤ n1/2║x║2 ≤ n║x║∞ 另外,若p和q是赫德尔(hölder)共轭指标,即1/p+1/q=1,那么有赫德尔不等式: | 一般来讲矩阵范数除了正定性,齐次性和三角不等式之外,还规定其必须满足相容性:║xy║≤║x║║y║。所以矩阵范数通常也称为相容范数。 如果║·║α是相容范数,且任何满足║·║β≤║·║α的范数║·║β都不是相容范数,那么║·║α称为极小范数。对于n阶实方阵(或复方阵)全体上的任何一个范数║·║,总存在唯一的实数k>0,使得k║·║是极小范数。 注:如果不考虑相容性,那么矩阵范数和向量范数就没有区别,因为mxn矩阵全体和mn维向量空间同构。引入相容性主要是为了保持矩阵作为线性算子的特征,这一点和算子范数的相容性一致,并且可以得到mincowski定理以外的信息。 第四章 算法总结 本学期讲解过的主要算法列举如下:线性方程组的解法(高斯消元法,列主消元法,doolittle分解法,追赶法,ldl分解法,jacobi分解法,seidel迭代法);非线性方程的求根方法(二分法,简单迭代法,newton迭代法,newton+下山因子,newton迭代法2,newton非线性方程);矩阵特征值与特征向量的计算(householder矩阵,反幂法,幂法,qr分解);函数的插值方法(三次样条插值,lagrange插值法,newton差商插值法);最佳平方逼近(chebyshev最小二乘法,曲线拟合最小二乘法);数值积分与数值微分(simpson求积分式算法,romberg算法,外推法);常微分方程初值问题的数值解法(欧拉改进法、龙格库塔法和修正的adams法)。下面对主要算法进行分析。4.1线性方程组的解法 本章学习了一些求解线性方程组的常用方法,其中gauss消元法,列主元消元法,lu分解法,追赶法和ldl’分解法都是解线性方程组的直接方法;而jacobi迭代法和sor法则是解线性方程组的基本迭代法。求解线性方程组时,应该注意方程组的性态,对病态方程组使用通常求解方程组的方法将导致错误。迭代求精法可用于求解某些病态方程。4.1.1高斯列主元lu分解法求解线性方程组 高斯消元法和lu分解法是直接法求解线性方程组中的两种方法。其中高斯消元法的基本思想是将线性方程组(1.1)通过消元,逐步化为同解的三角形方程组,然后用回代法解出n个解。高斯列主元消元法则是在高斯消元法的基础上提(k?1)(k?1)a?0akkkk出的先选主元再消元的方法,避免了时消元无法进行或者是当的绝(k?1)a(i?k?1,k?2,ik对值与其下方的元素,n)的绝对值之比很小时,引起计算机 上溢或产生很大的舍入误差而导致所求出的解失真的问题。lu分解法是将矩阵a用一个下三角矩阵和一个上三角矩阵之积来表示,即a?lu,然后由a?lu,ax?b,得lux?b,将线性方程组的求解化为对两个三角形方程组ly?b和ux?y的求解,由此可解出线性方程组(1.1)的n个解x1,x2,xn。这两种求解线性方程组的方法在处理单个线性方程组时没有差别,只是方法的不同,但在处理系数矩阵a相同,而右端项不同的一组线性方程组时,lu分解法就有明显的优势,因为它是将系数矩阵a和右端项b分开处理的,这样就可以只进行一次分解。例如,求解线性方程组ax?bi,i?1,2,m,用高斯消元法求解的计算量 1313mnn?mn2 大约为3,而用lu分解求解的计算量约为3,后者计算量显然小很多。但是lu分解法同样有可能由于ujj的绝对值很小而引起计算机上溢或产生很大的舍入误差而导致所求出的解失真。因此提出了结合高斯列主元消元的lu分解法。 我们采用的计算方法是先将a矩阵进行高斯列主元消元,然后再计算相应的l矩阵和u矩阵(u矩阵就是经过n-1步消元后的a矩阵)。但要注意,第k步消元时会产生mik(i?k?1,k?2,n),从而可以得到l矩阵的第k列元素,但在下一步消元前选取列主元时可能会交换方程的位置,因此与方程位置对应的l矩阵中的元素也要交换位置。4.2非线性方程组的求根方法 本章学习的二分法简单迭代法、newton迭代法等方法,代表着求解非线性方程所采用的两类方法。大范围收敛方法的初值x0选取没有多少限制,只要在含根区间任选其一即可,二分法就是这类方法。局部收敛法要求x0要充分靠近根x*才能保证收敛,以简单迭代法为基础,newton迭代法为代表的各类迭代法都属这类方法。4.2.1newton迭代法 牛顿迭代法的构造过程是这样的:设x0是f(x)?0的一个近似根,将f(x)在 f(x0)f(x)?f(x0)?f(x0)(x?x0)?(x?x0)2?x0处作taylor展开得2!,若取其 x?x?f(x)/f(x0),然后再对x1做f(x)100前两项来近似代替,得近似方程的根 f上述同样处理,继续下去,一般若(xk)?0,则可以构造出迭代格式 xk?1?xk?f(xk)f(xk)此格式称为牛顿迭代格式,用它来求解f(x)?0的方法称为牛顿迭代法。牛顿迭代法的几何意义是用f(x)在xk处的切线与x轴得交点作为下一个迭代点xk?1的。由于这一特点,牛顿迭代法也常称为切线法。 牛顿迭代法虽然收敛很快,但它通常过于依赖初值x0的选取,如果x0选择不当,将导致迭代发散或产生无限循环。 数值分析实验报告 一、实验3.1 题目: 考虑线性方程组,,编制一个能自动选取主元,又能手动选取主元的求解线性代数方程组的Gauss消去过程。 (1)取矩阵,则方程有解。取计算矩阵的条件数。分别用顺序Gauss消元、列主元Gauss消元和完全选主元Gauss消元方法求解,结果如何? (2)现选择程序中手动选取主元的功能,每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元,观察并记录计算结果,若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。 (3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。 (4)选取其他你感兴趣的问题或者随机生成的矩阵,计算其条件数,重复上述实验,观察记录并分析实验的结果。 1.算法介绍 首先,分析各种算法消去过程的计算公式,顺序高斯消去法: 第k步消去中,设增广矩阵中的元素(若等于零则可以判定系数矩阵为奇异矩阵,停止计算),则对k行以下各行计算,分别用乘以增广矩阵的第行并加到第行,则可将增广矩阵中第列中以下的元素消为零;重复此方法,从第1步进行到第n-1步,则可以得到最终的增广矩阵,即; 列主元高斯消去法: 第k步消去中,在增广矩阵中的子方阵中,选取使得,当时,对中第行与第行交换,然后按照和顺序消去法相同的步骤进行。重复此方法,从第1步进行第n-1步,就可以得到最终的增广矩阵,即; 完全主元高斯消去法: 第k步消去中,在增广矩阵中对应的子方阵中,选取使得,若或,则对中第行与第行、第列与第列交换,然后按照和顺序消去法相同的步骤进行即可。重复此方法,从第1步进行到第n-1步,就可以得到最终的增广矩阵,即; 接下来,分析回代过程求解的公式,容易看出,对上述任一种消元法,均有以下计算公式: 2.实验程序的设计 一、输入实验要求及初始条件; 二、计算系数矩阵A的条件数及方程组的理论解; 三、对各不同方法编程计算,并输出最终计算结果。 3.计算结果及分析 (1) 先计算系数矩阵的条件数,结果如下,可知系数矩阵的条件数较大,故此问题属于病态问题,b或A的扰动都可能引起解的较大误差; 采用顺序高斯消去法,计算结果为: 最终解为x=(1.***,1.***,1.***,1.***,0.***,1.***,0.***,1.***,0.***,1.***)T 使用无穷范数衡量误差,得到=2.842***1e-14,可以发现,采用顺序高斯消元法求得的解与精确解之间误差较小。通过进一步观察,可以发现,按照顺序高斯消去法计算时,其选取的主元值和矩阵中其他元素大小相近,因此顺序高斯消去法方式并没有对结果造成特别大的影响。 若采用列主元高斯消元法,则结果为: 最终解为x=(1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***)T 同样使用无穷范数衡量误差,有=0; 若使用完全主元高斯消元法,则结果为 最终解x=(1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***)T 同样使用无穷范数衡量误差,有=0; (2) 若每步都选取模最小或尽可能小的元素为主元,则计算结果为 最终解x=(1.*** 1.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.***)T 使用无穷范数衡量误差,有为2.842***1e-14;而完全主元消去法的误差为=0。 从(1)和(2)的实验结果可以发现,列主元消去法和完全主元消去法都得到了精确解,而顺序高斯消去法和以模尽量小的元素为主元的消去法没有得到精确解。在后两种消去法中,由于程序计算时的舍入误差,对最终结果产生了一定的影响,但由于方程组的维度较低,并且元素之间相差不大,所以误差仍比较小。 为进一步分析,计算上述4种方法每步选取的主元数值,并列表进行比较,结果如下: 第n次消元 顺序 列主元 完全主元 模最小 6.*** 6.*** 4.*** 4.*** 4.*** 4.*** 4.***3333 4.***3333 4.*** 4.*** 4.*** 4.*** 4.0***063 4.0***063 4.*** 4.*** 4.0039*** 4.0039*** 4.*** 0.0***469 0.0***469 4.*** 从上表可以发现,对这个方程组而言,顺序高斯消去选取的主元恰好事模尽量小的元素,而由于列主元和完全主元选取的元素为8,与4在数量级上差别小,所以计算过程中的累积误差也较小,最终4种方法的输出结果均较为精确。 在这里,具体解释一下顺序法与模最小法的计算结果完全一致的原因。该矩阵在消元过程中,每次选取主元的一列只有两个非零元素,对角线上的元素为4左右,而其正下方的元素为8,该列其余位置的元素均为0。在这样的情况下,默认的主元也就是该列最小的主元,因此两种方法所得到的计算结果是一致的。 理论上说,完全高斯消去法的误差最小,其次是列主元高斯消去法,而选取模最小的元素作为主元时的误差最大,但是由于方程组的特殊性(元素相差不大并且维度不高),这个理论现象在这里并没有充分体现出来。 (3) 时,重复上述实验过程,各种方法的计算结果如下所示,在这里,仍采用无穷范数衡量绝对误差。 顺序高斯消去法 列主元高斯消去 完全主元高斯消去 选取模最小或尽可能小元素作为主元消去 X 1.*** 1.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 1.*** 1.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 2.***e-11 0 0 2.***e-11 可以看出,此时列主元和完全主元的计算结果仍为精确值,而顺序高斯消去和模尽可能小方法仍然产生了一定的误差,并且两者的误差一致。与n=10时候的误差比相比,n=20时的误差增长了大约1000倍,这是由于计算过程中舍入误差的不断累积所致。所以,如果进一步增加矩阵的维数,应该可以看出更明显的现象。 (4) 不同矩阵维度下的误差如下,在这里,为方便起见,选取2-条件数对不同维度的系数矩阵进行比较。 维度 条件数 顺序消去 列主元 完全主元 模尽量小 1.7e+3 2.84e-14 0 0 2.84e-14 1.8e+6 2.91e-11 0 0 2.91e-11 5.7e+7 9.31e-10 0 0 9.31e-10 1.8e+9 2.98e-08 0 0 2.98e-08 1.9e+12 3.05e-05 0 0 3.05e-05 3.8e+16 3.28e+04 3.88e-12 3.88e-12 3.28e+04 8.5e+16 3.52e+13 4.2e-3 4.2e-3 3.52e+13 从上表可以看出,随着维度的增加,不同方法对计算误差的影响逐渐体现,并且增长较快,这是由于舍入误差逐步累计而造成的。不过,方法二与方法三在维度小于40的情况下都得到了精确解,这两种方法的累计误差远比方法一和方法四慢;同样地,出于与前面相同的原因,方法一与方法四的计算结果保持一致,方法二与方法三的计算结果保持一致。 4.结论 本文矩阵中的元素差别不大,模最大和模最小的元素并没有数量级上的差异,因此,不同的主元选取方式对计算结果的影响在维度较低的情况下并不明显,四种方法都足够精确。 对比四种方法,可以发现采用列主元高斯消去或者完全主元高斯消去法,可以尽量抑制误差,算法最为精确。不过,对于低阶的矩阵来说,四种方法求解出来的结果误差均较小。 另外,由于完全选主元方法在选主元的过程中计算量较大,而且可以发现列主元法已经可以达到很高的精确程度,因而在实际计算中可以选用列主元法进行计算。 附录:程序代码 clear clc; format long; %方法选择 n=input('矩阵A阶数:n='); disp('选取求解方式'); disp('1 顺序Gauss消元法,2 列主元Gauss消元法,3 完全选主元Gauss消元法,4 模最小或近可能小的元素作为主元'); a=input('求解方式序号:'); %赋值A和b A=zeros(n,n); b=zeros(n,1); for i=1:n A(i,i)=6; if i>1 A(i,i-1)=8; end if i A(i,i+1)=1; end end for i=1:n for j=1:n b(i)=b(i)+A(i,j); end end disp('给定系数矩阵为:'); A disp('右端向量为:'); b %求条件数及理论解 disp('线性方程组的精确解:'); X=(A\b)' fprintf('矩阵A的1-条件数: %f \n',cond(A,1)); fprintf('矩阵A的2-条件数: %f \n',cond(A)); fprintf('矩阵A的无穷-条件数: %f \n',cond(A,inf)); %顺序Gauss消元法 if a==1 A1=A;b1=b; for k=1:n if A1(k,k)==0 disp('主元为零,顺序Gauss消元法无法进行'); break end fprintf('第%d次消元所选取的主元:%g\n',k,A1(k,k)) %disp('此次消元后系数矩阵为:'); %A1 for p=k+1:n l=A1(p,k)/A1(k,k); A1(p,k:n)=A1(p,k:n)-l*A1(k,k:n); b1(p)=b1(p)-l*b1(k); end end x1(n)=b1(n)/A1(n,n); for k=n-1:-1:1 for w=k+1:n b1(k)=b1(k)-A1(k,w)*x1(w); end x1(k)=b1(k)/A1(k,k); end disp('顺序Gauss消元法解为:'); disp(x1); disp('所求解与精确解之差的无穷-范数为'); norm(x1-X,inf) end %列主元Gauss消元法 if a==2 A2=A;b2=b; for k=1:n [max_i,max_j]=find(A2(:,k)==max(abs(A2(k:n,k)))); if max_i~=k A2_change=A2(k,:); A2(k,:)=A2(max_i,:); A2(max_i,:)=A2_change; b2_change=b2(k); b2(k)=b2(max_i); b2(max_i)=b2_change; end if A2(k,k)==0 disp('主元为零,列主元Gauss消元法无法进行'); break end fprintf('第%d次消元所选取的主元:%g\n',k,A2(k,k)) %disp('此次消元后系数矩阵为:'); %A2 for p=k+1:n l=A2(p,k)/A2(k,k); A2(p,k:n)=A2(p,k:n)-l*A2(k,k:n); b2(p)=b2(p)-l*b2(k); end end x2(n)=b2(n)/A2(n,n); for k=n-1:-1:1 for w=k+1:n b2(k)=b2(k)-A2(k,w)*x2(w); end x2(k)=b2(k)/A2(k,k); end disp('列主元Gauss消元法解为:'); disp(x2); disp('所求解与精确解之差的无穷-范数为'); norm(x2-X,inf) end %完全选主元Gauss消元法 if a==3 A3=A;b3=b; for k=1:n VV=eye(n); [max_i,max_j]=find(A3(k:n,k:n)==max(max(abs(A3(k:n,k:n))))); if numel(max_i)==0 [max_i,max_j]=find(A3(k:n,k:n)==-max(max(abs(A3(k:n,k:n))))); end W=eye(n); W(max_i(1)+k-1,max_i(1)+k-1)=0; W(k,k)=0; W(max_i(1)+k-1,k)=1; W(k,max_i(1)+k-1)=1; V=eye(n); V(k,k)=0; V(max_j(1)+k-1,max_j(1)+k-1)=0; V(k,max_j(1)+k-1)=1; V(max_j(1)+k-1,k)=1; A3=W*A3*V; b3=W*b3; VV=VV*V; if A3(k,k)==0 disp('主元为零,完全选主元Gauss消元法无法进行'); break end fprintf('第%d次消元所选取的主元:%g\n',k,A3(k,k)) %disp('此次消元后系数矩阵为:'); %A3 for p=k+1:n l=A3(p,k)/A3(k,k); A3(p,k:n)=A3(p,k:n)-l*A3(k,k:n); b3(p)=b3(p)-l*b3(k); end end x3(n)=b3(n)/A3(n,n); for k=n-1:-1:1 for w=k+1:n b3(k)=b3(k)-A3(k,w)*x3(w); end x3(k)=b3(k)/A3(k,k); end disp('完全选主元Gauss消元法解为:'); disp(x3); disp('所求解与精确解之差的无穷-范数为'); norm(x3-X,inf) end %模最小或近可能小的元素作为主元 if a==4 A4=A;b4=b; for k=1:n AA=A4; AA(AA==0)=NaN; [min_i,j]=find(AA(k:n,k)==min(abs(AA(k:n,k)))); if numel(min_i)==0 [min_i,j]=find(AA(k:n,k)==-min(abs(AA(k:n,k:n)))); end W=eye(n); W(min_i(1)+k-1,min_i(1)+k-1)=0; W(k,k)=0; W(min_i(1)+k-1,k)=1; W(k,min_i(1)+k-1)=1; A4=W*A4; b4=W*b4; if A4(k,k)==0 disp('主元为零,模最小Gauss消元法无法进行'); break end fprintf('第%d次消元所选取的主元:%g\n',k,A4(k,k)) %A4 for p=k+1:n l=A4(p,k)/A4(k,k); A4(p,k:n)=A4(p,k:n)-l*A4(k,k:n); b4(p)=b4(p)-l*b4(k); end end x4(n)=b4(n)/A4(n,n); for k=n-1:-1:1 for w=k+1:n b4(k)=b4(k)-A4(k,w)*x4(w); end x4(k)=b4(k)/A4(k,k); end disp('模最小Gauss消元法解为:'); disp(x4); disp('所求解与精确解之差的无穷-范数为'); norm(x4-X,inf) end 二、实验3.3 题目: 考虑方程组的解,其中系数矩阵H为Hilbert矩阵: 这是一个著名的病态问题。通过首先给定解(例如取为各个分量均为1)再计算出右端的办法给出确定的问题。 (1)选择问题的维数为6,分别用Gauss消去法(即LU分解)、J迭代法、GS迭代法和SOR迭代法求解方程组,其各自的结果如何?将计算结果与问题的解比较,结论如何。 (2)逐步增大问题的维数,仍用上述的方法来解它们,计算的结果如何?计算的结果说明的什么? (3)讨论病态问题求解的算法。 1.算法设计 对任意线性方程组,分析各种方法的计算公式如下,(1)Gauss消去法: 首先对系数矩阵进行LU分解,有,则原方程转化为,令,则原方程可以分为两步回代求解: 具体方法这里不再赘述。 (2)J迭代法: 首先分解,再构造迭代矩阵,其中,进行迭代计算,直到误差满足要求。 (3)GS迭代法: 首先分解,再构造迭代矩阵,其中,进行迭代计算,直到误差满足要求。 (4)SOR迭代法: 首先分解,再构造迭代矩阵,其中,进行迭代计算,直到误差满足要求。 2.实验过程 一、根据维度n确定矩阵H的各个元素和b的各个分量值; 二、选择计算方法(Gauss消去法,J迭代法,GS迭代法,SOR迭代法),对迭代法设定初值,此外SOR方法还需要设定松弛因子; 三、进行计算,直至满足误差要求(对迭代法,设定相邻两次迭代结果之差的无穷范数小于0.0001; 3.计算结果及分析 (1)时,问题可以具体定义为 计算结果如下,Gauss消去法 第1次消元所选取的主元是:1 第2次消元所选取的主元是:0.0833333 第3次消元所选取的主元是:0.00555556 第4次消元所选取的主元是:0.000357143 第5次消元所选取的主元是:2.26757e-05 第6次消元所选取的主元是:1.43155e-06 解得X=(0.*** 1.*** 0.*** 1.*** 0.*** 1.***)T 使用无穷范数衡量误差,可得=4.***e-10; J迭代法 设定迭代初值为零,计算得到 J法的迭代矩阵B的谱半径为4.30853>1,所以J法不收敛; GS迭代法 设定迭代初值为零,计算得到GS法的迭代矩阵G的谱半径为:0.999998<1,故GS法收敛,经过541次迭代计算后,结果为X=(1.001***6 0.*** 0.*** 1.*** 1.*** 0.***)T 使用无穷范数衡量误差,有=0.***; SOR迭代法 设定迭代初值为零向量,并设定,计算得到SOR法迭代矩阵谱半径为0.***,经过100次迭代后的计算结果为 X=(1.*** 0.*** 1.03***59 1.06***81 1.*** 0.9***527)T; 使用无穷范数衡量误差,有=0.***; 对SOR方法,可变,改变值,计算结果可以列表如下 迭代次数 迭代矩阵的谱半径 0.*** 0.*** 0.*** 0.*** X 1.*** 0.*** 1.01***40 1.*** 1.0***681 0.*** 1.*** 0.*** 1.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.*** 1.*** 0.*** 0.*** 1.05***66 0.*** 1.*** 0.*** 1.*** 0.*** 0.*** 0.*** 0.*** 0.*** 可以发现,松弛因子的取值对迭代速度造成了不同的影响,上述四种方法中,松弛因子=0.5时,收敛相对较快。 综上,四种算法的结果列表如下: 算法 Gauss消去法 Jacobi法 GS法 SOR法(取) 迭代次数 -- 不收敛 541 迭代矩阵的谱半径 -- 4.30853 0.999998 0.*** X 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** -- 1.001***6 0.*** 0.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.03***59 1.06***81 1.*** 0.9***527 4.***e-10 -- 0.*** 0.*** 计算可得,矩阵H的条件数为>>1,所以这是一个病态问题。由上表可以看出,四种方法的求解都存在一定的误差。下面分析误差的来源: LU分解方法的误差存在主要是由于Hilbert矩阵各元素由分数形式转换为小数形式时,不能除尽情况下会出现舍入误差,在进行LU分解时也存在这个问题,所以最后得到的结果不是方程的精确解,但结果显示该方法的误差非常小; Jacobi迭代矩阵的谱半径为4.30853,故此迭代法不收敛; GS迭代法在迭代次数为541次时得到了方程的近似解,其误差约为0.05,比较大。GS迭代矩阵的谱半径为0.999998,很接近1,所以GS迭代法收敛速度较慢; SOR迭代法在迭代次数为100次时误差约为0.08,误差较大。SOR迭代矩阵的谱半径为0.999999,也很接近1,所以时SOR迭代法收敛速度不是很快,但是相比于GS法,在迭代速度方面已经有了明显的提高;另外,对不同的,SOR方法的迭代速度会相应有变化,如果选用最佳松弛因子,可以实现更快的收敛; (2) 考虑不同维度的情况,时,算法 Gauss消去 J法 GS法 SOR法(w=0.5) 计算结果 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** -- 0.*** 1.*** 0.*** 1.*** 1.*** 1.*** 0.9968*** 0.*** 1.*** 0.9397*** 0.*** 1.*** 1.*** 1.*** 0.*** 0.*** 迭代次数 -- -- 356 谱半径 -- 6.04213 0.*** -- 时,算法 Gauss消去法 Jacobi法 GS法 SOR法(w=0.5) 计算结果 0.*** 1.*** 0.*** 1.000***1 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** -- 0.*** 1.*** 0.*** 0.*** 0.*** 1.02***91 1.*** 1.*** 1.*** 0.*** 0.947***7 1.0***572 0.*** 0.*** 0.*** 1.*** 1.*** 1.*** 1.*** 0.*** 0.*** 0.*** 迭代次数 -- -- 1019 谱半径 -- 8.64964 0.*** -- 时 算法 Gauss消去法 Jacobi法 GS法 SOR法(w=0.5) 计算结果 0.*** 1.*** 0.*** 0.*** 1.*** 0.*** 2.*** -2.*** 7.*** -7.*** 7.*** -1.*** 0.*** 1.*** 0.*** -- 不收敛 1.*** 1.*** 0.907***9 0.*** 0.*** 1.*** 1.09***64 1.*** 1.*** 1.*** 1.0385*** 0.*** 0.942***3 0.*** 0.*** 迭代次数 -- -- 262 谱半径 -- 6.04213 >1 1.*** 8.*** -- -- 0.*** 分析以上结果可以发现,随着n值的增加,Gauss消去法误差逐渐增大,而且误差增大的速度很快,在维数小于等于10情况下,Gauss消去法得到的结果误差较小;但当维数达到15时,计算结果误差已经达到精确解的很多倍; J法迭代不收敛,无论n如何取值,其谱半径始终大于1,因而J法不收敛,所以J迭代法不能用于Hilbert矩阵的求解; 对于GS迭代法和SOR迭代法,两种方法均收敛,GS迭代法是SOR迭代法松弛因子取值为1的特例,SOR方法受到取值的影响,会有不同的收敛情况。可以得出GS迭代矩阵的谱半径小于1但是很接近1,收敛速度很慢。虽然随着维数的增大,所需迭代的次数逐渐减少,但是当维数达到15的时候,GS法已经不再收敛。因此可以得出结论,GS迭代方法在Hilbert矩阵维数较低时,能够在一定程度上满足迭代求解的需求,不过迭代的速度很慢。另外,随着矩阵维数的增加,SOR法的误差水平基本稳定,而且误差在可以接受的范围之内。 经过比较可以得出结论,如果求解较低维度的Hibert矩阵问题,Gauss消去法、GS迭代法和SOR迭代法均可使用,且Gauss消去法的结果精确度较高;如果需要求解较高维度的Hibert矩阵问题,只有采用SOR迭代法。 (3) 系数矩阵的条件数较大时,为病态方程。由实验可知,Gauss法在解上述方程时,结果存在很大的误差。而对于收敛的迭代法,可以通过选取最优松弛因子的方法来求解,虽然迭代次数相对较多,但是结果较为精确。 总体来看,对于一般病态方程组的求解,可以采用以下方式: 1.低维度下采用Gauss消去法直接求解是可行的; Jacobi迭代方法不适宜于求解病态问题; GS迭代方法可以解决维数较低的病态问题,但其谱半径非常趋近于1,导致迭代算法收敛速度很慢,维数较大的时候,GS法也不再收敛; SOR方法较适合于求解病态问题,特别是矩阵维数较高的时候,其优势更为明显。 2.采用高精度的运算,如选用双倍或更多倍字长的运算,可以提高收敛速度; 3.可以对原方程组作某些预处理,从而有效降低系数矩阵的条件数。 4.实验结论 (1)对Hibert矩阵问题,其条件数会随着维度的增加迅速增加,病态性会越来越明显;在维度较低的时候,Gauss消去法、GS迭代法和SOR迭代法均可使用,且可以优先使用Gauss消去法;如果需要求解较高维度的Hibert矩阵问题,只有SOR迭代法能够求解。 (2)SOR方法比较适合于求解病态问题,特别是矩阵维数较高的时候,其优点更为明显。从本次实验可以看出,随着矩阵维数的增大,SOR方法所需的迭代次数减少,而且误差基本稳定,是解决病态问题的适宜方法。 附录:程序代码 clear all clc; format long; %矩阵赋值 n=input('矩阵H的阶数:n='); for i=1:n for j=1:n H(i,j)=1/(i+j-1); end end b=H*ones(n,1); disp('H矩阵为:'); H disp('向量b:'); b %方法选择 disp('选取求解方式'); disp('1 Gauss消去法,2 J迭代法,3 GS迭代法,4 SOR迭代法'); a=input('求解方式序号:'); %Gauss消去法 if a==1; H1=H;b1=b; for k=1:n if H1(k,k)==0 disp('主元为零,Gauss消去法无法进行'); break end fprintf('第%d次消元所选取的主元是:%g\n',k,H1(k,k)) for p=k+1:n m5=-H1(p,k)/H1(k,k); H1(p,k:n)=H1(p,k:n)+m5*H1(k,k:n); b1(p)=b1(p)+m5*b1(k); end end x1(n)=b1(n)/H1(n,n); for k=n-1:-1:1 for v=k+1:n b1(k)=b1(k)-H1(k,v)*x1(v); end x1(k)=b1(k)/H1(k,k); end disp('Gauss消去法解为:'); disp(x1); disp('解与精确解之差的无穷范数'); norm((x1-a),inf) end D=diag(diag(H)); L=-tril(H,-1); U=-triu(H,1); %J迭代法 if a==2; %给定初始x0 ini=input('初始值设定:x0='); x0(:,1)=ini*diag(ones(n)); disp('初始解向量为:'); x0 xj(:,1)=x0(:,1); B=(D^(-1))*(L+U); f=(D^(-1))*b; fprintf('(J法B矩阵谱半径为:%g\n',vrho(B)); if vrho(B)<1; for m2=1:5000 xj(:,m2+1)=B*xj(:,m2)+fj; if norm((xj(:,m2+1)-xj(:,m2)),inf)<0.0001 break end end disp('J法计算结果为:'); xj(:,m2+1) disp('解与精确解之差的无穷范数'); norm((xj(:,m2+1)-diag(ones(n))),inf) disp('J迭代法迭代次数:'); m2 else disp('由于B矩阵谱半径大于1,因而J法不收敛'); end end %GS迭代法 if a==3; %给定初始x0 ini=input('初始值设定:x0='); x0(:,1)=ini*diag(ones(n)); disp('初始解向量为:'); x0 xG(:,1)=x0(:,1); G=inv(D-L)*U; fG=inv(D-L)*b; fprintf('GS法G矩阵谱半径为:%g\n',vrho(G)); if vrho(G)<1 for m3=1:5000 xG(:,m3+1)=G*xG(:,m3)+fG; if norm((xG(:,m3+1)-xG(:,m3)),inf)<0.0001 break; end end disp('GS迭代法计算结果:'); xG(:,m3+1) disp('解与精确解之差的无穷范数'); norm((xG(:,m3+1)-diag(ones(n))),inf) disp('GS迭代法迭代次数:'); m3 else disp('由于G矩阵谱半径大于1,因而GS法不收敛'); end end %SOR迭代法 if a==4; %给定初始x0 ini=input('初始值设定:x0='); x0(:,1)=ini*diag(ones(n)); disp('初始解向量为:'); x0 A=H; for i=1:n b(i)=sum(A(i,:)); end x_star=ones(n,1); format long w=input('松弛因子:w='); Lw=inv(D-w*L)*((1-w)*D+w*U); f=w*inv(D-w*L)*b; disp('迭代矩阵的谱半径:') p=vrho(Lw) time_max=100;%迭代次数 x=zeros(n,1);%迭代初值 for i=1:time_max x=Lw*x+f; end disp('SOR迭代法得到的解为'); x disp('解与精确解之差的无穷范数'); norm((x_star-x),inf) end pause 三、实验4.1 题目: 对牛顿法和拟牛顿法。进行非线性方程组的数值求解 (1)用上述两种方法,分别计算下面的两个例子。在达到精度相同的前提下,比较其迭代次数、CPU时间等。 (2)取其他初值,结果又如何?反复选取不同的初值,比较其结果。 (3)总结归纳你的实验结果,试说明各种方法适用的问题。 1.算法设计 对需要求解的非线性方程组而言,牛顿法和拟牛顿法的迭代公式如下,(1)牛顿法: 牛顿法为单步迭代法,需要取一个初值。 (2)拟牛顿法:(Broyden秩1法) 其中,拟牛顿法不需要求解的导数,因此节省了大量的运算时间,但需要给定矩阵的初值,取为。 2.实验过程 一、输入初值; 二、根据误差要求,按公式进行迭代计算; 三、输出数据; 3.计算结果及分析 (1)首先求解方程组(1),在这里,设定精度要求为,方法 牛顿法 拟牛顿法 初始值 计算结果X x1 0.*** 0.*** x2 1.*** 1.0852*** x3 0.*** 0.*** 迭代次数 CPU计算时间/s 3.777815 2.739349 可以看出,在初始值相同情况下,牛顿法和拟牛顿法在达到同样计算精度情况下得到的结果基本相同,但牛顿法的迭代次数明显要少一些,但是,由于每次迭代都需要求解矩阵的逆,所以牛顿法每次迭代的CPU计算时间更长。 之后求解方程组(2),同样设定精度要求为 方法 牛顿法 拟牛顿法 初始值 计算结果X x1 0.*** 0.*** x2 0.*** 0.*** x3 -0.*** -0.*** 迭代次数 CPU计算时间/s 2.722437 3.920195 同样地,可以看出,在初始值相同情况下,牛顿法和拟牛顿法在达到同样计算精度情况下得到的结果是基本相同的,但牛顿法的迭代次数明显要少,但同样的,由于每次迭代中有求解矩阵的逆的运算,牛顿法每次迭代的CPU计算时间较长。 (2)对方程组(1),取其他初值,计算结果列表如下,同样设定精度要求为 初始值 方法 牛顿法 拟牛顿法 计算结果 0.*** 1.*** 0.*** 9.21***94 -5.*** 18.1***205 迭代次数 CPU计算时间/s 3.907164 4.818019 计算结果 0.*** 1.*** 0.*** 9.21***91 -5.*** 18.1***807 迭代次数 2735 CPU计算时间/s 8.127286 5.626023 计算结果 0.*** 1.*** 0.*** 0.*** 1.0852*** 0.*** 迭代次数 CPU计算时间/s 3.777815 2.739349 计算结果 0.*** 1.*** 0.*** 0.*** 1.*** 0.*** 迭代次数 188 CPU计算时间/s 3.835697 2.879070 计算结果 9.21***22 -5.*** 18.1***605 Matlab警告矩阵接近奇异值,程序进入长期循环计算中 迭代次数 -- CPU计算时间/s 4.033868 -- 计算结果 0.*** 1.*** 0.*** Matlab警告矩阵接近奇异值,程序进入长期循环计算中 迭代次数 -- CPU计算时间/s 12.243263 -- 从上表可以发现,方程组(1)存在另一个在(9.2,-5.6,18.1)T附近的不动点,初值的选取会直接影响到牛顿法和拟牛顿法最后的收敛点。 总的来说,设定的初值离不动点越远,需要的迭代次数越多,因而初始值的选取非常重要,合适的初值可以更快地收敛,如果初始值偏离精确解较远,会出现迭代次数增加直至无法收敛的情况; 由于拟牛顿法是一种近似方法,拟牛顿法需要的的迭代次数明显更多,而且收敛情况不如牛顿法好(初值不够接近时,甚至会出现奇异矩阵的情况),但由于牛顿法的求解比较复杂,计算时间较长; 同样的,对方程组(2),取其他初值,计算结果列表如下,同样设定精度要求为 初始值 方法 牛顿法 拟牛顿法 计算结果 0.*** 0.*** -0.*** 0.*** 0.*** -0.*** 迭代次数 CPU计算时间/s 2.722437 3.920195 计算结果 0.*** 0.*** -0.*** 0.*** -0.*** 76.*** 迭代次数 CPU计算时间/s 5.047111 5.619752 计算结果 0.*** 0.*** -0.*** 1.0e+02 * -0.*** -0.000***6 1.754***3 迭代次数 CPU计算时间/s 3.540668 3.387829 计算结果 0.*** 0.*** -0.*** 1.0e+04 * 0.*** -0.*** 1.*** 迭代次数 CPU计算时间/s 2.200571 2.640901 计算结果 0.*** 0.*** -0.*** 矩阵为奇异值,无法输出准确结果 迭代次数 -- CPU计算时间/s 1.719072 -- 计算结果 0.*** 0.*** -0.*** 矩阵为奇异值,无法输出准确结果 迭代次数 149 -- CPU计算时间/s 2.797116 -- 计算结果 矩阵为奇异值,无法输出准确结果 矩阵为奇异值,无法输出准确结果 迭代次数 -- -- CPU计算时间/s -- -- 在这里,与前文类似的发现不再赘述。 从这里看出,牛顿法可以在更大的区间上实现压缩映射原理,可以在更大的范围上选取初值并最终收敛到精确解附近; 在初始值较接近于不动点时,牛顿法和拟牛顿法计算所得到的结果是基本相同的,虽然迭代次数有所差别,但计算总的所需时间相近。 (3) 牛顿法在迭代过程中用到了矩阵的求逆,其迭代收敛的充分条件是迭代满足区间上的映内性,对于矩阵的求逆过程比较简单,所以在较大区间内满足映内性的问题适合应用牛顿法进行计算。一般而言,对于函数单调或者具有单值特性的函数适合应用牛顿法,其对初始值敏感程度较低,算法具有很好的收敛性。 另外,需要说明的是,每次计算给出的CPU时间与计算机当时的运行状态有关,同时,不同代码的运行时间也不一定一致,所以这个数据并不具有很大的参考价值。 4.实验结论 对牛顿法和拟牛顿法,都存在初始值越接近精确解,所需的迭代次数越小的现象; 在应用上,牛顿法和拟牛顿法各有优势。就迭代次数来说,牛顿法由于更加精确,所需的迭代次数更少;但就单次迭代来说,牛顿法由于计算步骤更多,且计算更加复杂,因而每次迭代所需的时间更长,而拟牛顿法由于采用了简化的近似公式,其每次迭代更加迅速。当非线性方程组求逆过程比较简单时,如方程组1的情况时,拟牛顿法不具有明显的优势;而当非线性方程组求逆过程比较复杂时,如方程组2的情况,拟牛顿法就可以体现出优势,虽然循环次数有所增加,但是CPU耗时反而更少。 另外,就方程组压缩映射区间来说,一般而言,对于在区间内函数呈现单调或者具有单值特性的函数适合应用牛顿法,其对初始值敏感程度较低,使算法具有很好的收敛性;而拟牛顿法由于不需要在迭代过程中对矩阵求逆,而是利用差商替代了对矩阵的求导,所以即使初始误差较大时,其倒数矩阵与差商偏差也较小,所以对初始值的敏感程度较小。 附录:程序代码 %方程1,牛顿法 tic; format long; %%初值 disp('请输入初值'); a=input('第1个分量为:'); b=input('第2个分量为:'); c=input('第3个分量为:'); disp('所选定初值为'); x=[a;b;c] %%误差要求 E=0.0001; %%迭代 i=0; e=2*E; while e>E F=[12*x(1)-x(2)^2-4*x(3)-7;x(1)^2+10*x(2)-x(3)-11;x(2)^3+10*x(3)-8]; f=[12,-2*x(2),-4;2*x(1),10,-1;0,3*x(2)^2,10]; det_x=((f)^(-1))*(-F); x=x+det_x; e=max(norm(det_x)); i=i+1; end disp('迭代次数'); i disp('迭代次数'); x toc; %方程1,拟牛顿法 tic; format long; %%初值 %%初值 disp('请输入初值'); a=input('第1个分量为:'); b=input('第2个分量为:'); c=input('第3个分量为:'); disp('所选定初值为'); x0=[a;b;c] %%误差要求 E=0.0001; %%迭代 i=0; e=2*E; A0=eye(3); while e>E F0=[12*x0(1)-x0(2)^2-4*x0(3)-7;x0(1)^2+10*x0(2)-x0(3)-11;x0(2)^3+10*x0(3)-8]; x1=x0-A0^(-1)*F0; s=x1-x0; F1=[12*x1(1)-x1(2)^2-4*x1(3)-7;x1(1)^2+10*x1(2)-x1(3)-11;x1(2)^3+10*x1(3)-8]; y=F1-F0; A1=A0+(y-A0*s)*s'/(s'*s); x0=x1; A0=A1; e=max(norm(s)); i=i+1; end disp('迭代次数'); i disp('迭代次数'); x0 toc; %方程2,牛顿法 tic; format long; %%初值 disp('请输入初值'); a=input('第1个分量为:'); b=input('第2个分量为:'); c=input('第3个分量为:'); disp('所选定初值为'); x=[a;b;c] %%误差要求 E=0.0001; %%迭代 i=0; e=2*E; while e>E F=[3*x(1)-cos(x(2)*x(3))-0.5;x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;exp(1)^(-x(1)*x(2))+20*x(3)+(10*pi-3)/3]; f=[3,x(3)*sin(x(2)*x(3)),x(2)*sin(x(2)*x(3));2*x(1),-162*x(2)-81/5,cos(x(3));-x(2)*exp(1)^(-x(1)*x(2)),-x(1)*exp(1)^(-x(1)*x(2)),20]; det_x=((f)^(-1))*(-F); x=x+det_x; e=max(norm(det_x)); i=i+1; end disp('迭代次数'); i disp('迭代次数'); x toc; %方程2,拟牛顿法 tic; format long; %%初值 %%初值 disp('请输入初值'); a=input('第1个分量为:'); b=input('第2个分量为:'); c=input('第3个分量为:'); disp('所选定初值为'); x0=[a;b;c] %%误差要求 E=0.0001; %%迭代 i=0; e=2*E; A0=eye(3); while e>E F0=[3*x0(1)-cos(x0(2)*x0(3))-0.5;x0(1)^2-81*(x0(2)+0.1)^2+sin(x0(3))+1.06;exp(1)^(-x0(1)*x0(2))+20*x0(3)+(10*pi-3)/3]; x1=x0-A0^(-1)*F0; s=x1-x0; F1=[3*x1(1)-cos(x1(2)*x1(3))-0.5;x1(1)^2-81*(x1(2)+0.1)^2+sin(x1(3))+1.06;exp(1)^(-x1(1)*x1(2))+20*x1(3)+(10*pi-3)/3]; y=F1-F0; A1=A0+(y-A0*s)*s'/(s'*s); x0=x1; A0=A1; e=max(norm(s)); i=i+1; end disp('迭代次数'); i disp('迭代次数'); x0 toc; 2011级葫芦岛校区研究生数值分析复习参考提纲(注意例题未必出原题,给出的是题型) 一、例2-4,例2-6,例2-11,二、86页:1,2,3,5,6,7,8 三、1、n阶线性方程组的雅可比迭代法:迭代公式、矩阵表示 2、n阶线性方程组的高斯-赛德尔迭代法:迭代公式、矩阵表示 3、逐次超松弛迭代法:迭代公式、矩阵表示 4、迭代法的收敛性:(1)迭代法收敛的充要条件;(2)迭代法收敛的充分条件;(3)迭代法收敛的基本定理 5、例4-3,例4-4,例4-5,例4-8,6、110页:1,2,5 四、1、n次拉格朗日插值:插值公式、插值余项与误差估计 2、牛顿插值:差商、牛顿插值公式、插值余项与误差估记 3、等距节点插值:差分、等距节点插值公式、误差估计 4、厄米特插值:插值公式、插值余项 5、分段低次插值:插值方法,分段低次线性插值余项、分段三次厄米特插值余项 6、三次样条插值:问题的提法、边界条件类型 7、例5- 1、例5- 2、例5- 4、例5- 6、例5- 7、例5- 8、例5-10 五、1、最佳一致逼近的概念、最佳一致逼近多项式的解法、2、切比雪夫多项式及其性质 3、例6-24、最佳平方逼近的概念、最佳平方逼近的计算 5、例6-56、正交多项式序列的构造方法 7、勒让德多项式及其性质 8、例6-69、离散数据的最小二乘法 10、例6-8 六、1、牛顿-柯特斯求积公式、梯形公式、辛普森公式、柯特斯公式及其阶段误差 2、证明:梯形公式和矩形公式具有一次代数精度、辛普森公式具有三次代数精度 3、181页定理1及其证明 4、复化梯形公式、复化辛普森公式,复化柯特斯公式及其截断误差 5、例7- 2、例7-36、龙贝格求积公式 7、例7- 5、例7-68、插值型求导方法 9、例7-12 七、例8- 1、例8-2 八、10页1、3、7、9 第六章 数值积分 --------学习小结 姓名 班级 学号 一、本章学习体会 本章主要讲授了数值积分的一些求积公式及各种求积公式的代数精度,重点应掌握插值型求积公式,什么样的求积公式可以被称为插值型求积公式,Newton-Cotes求积公式及其收敛性与数值稳定性,复化求积公式和高斯求积公式,在本章的学习过程中也遇到不少问题,比如本章知识点多,公式多,在做题时容易张冠李戴,其次对Newton-Cotes求积公式的收敛性与数值稳定性理解不够透彻,处理一个实际问题时,不知道选取哪一种求积公式,来达到最精确的结果。 二、本章知识梳理 6.1求积公式及其代数精度 代数精度的概念:如果求积公式(6.1)当f(x)为任何次数不高于m的多项式时都成为等式,而当f(x)为某个m+1次多项式时(6.1)不能成为等式,则称求积公式(6.1)具有m次代数精度。6.2插值型求积公式 (1)求积公式: Rnabf(n1)()n1(x)dx (n1)!(2)重要的定理:n+1个节点的插值型求积公式至少具有n次代数度。(3)求积系数: k0nAkba 6.3Newton-Cotes求积公式及其收敛性与数值稳定性 (n)f(xk)(1)公式:f(x)dxf(xk)(ba)cka(n)kk0k0bnnnhn2n(n1)(2)截断误差:Rnf()(ttj)dt (n1)!0j0(3)重要的定理:当n为偶数时,n+1个节点的Newton-Cotes求积公式至少具有n+1次代数精度。 (4)常用的Newton-Cotes求积公式 n=1 梯形公式:babaf(x)dx[f(a)f(b)] 2(ba)3f(),(a,b),具有一次精度。 余项:R112n=2 Simpson公式:f(x)dxabbaab[f(a)4f()f(b)] 62(ba)5(4)f(),(a,b),具有三次精度。余项:R228806.4复化求积法 (1)复化梯形公式: 截断误差: ban1hf(x)dx[f(a)f(b)2f(akh)]2k1 RTba2hf(),[a,b]12 (2)复化Simpson公式: bamm1hf(x)dx[f(a)f(b)4f(x2k1)2f(x2k)]3k1k1 截断误差: Rsba4(4)hf(),[a,b]180 6.5Gauss型求积公式 (1)定义:若n个节点的插值型求积公式(6.23)具有2n-1 次代数精度,则称它为Gauss型求积公式。 (2)定理:n个节点的 Gauss型求积公式的代数精度为2n-1。 (3)定理:设{gk(x),k0,1,}是区间[a,b]上带权(x)的正交多项式系,则求积公式(6.23)、式(6.24)是Gauss型求积公式的充分必要条件是它的求积节点是n次正交多项式gn(x)的n个零点。(4)求积系数 公式: Akb(x)gn(x)(xk)(xxk)gnadx,k1,2,,n 性质:1.Ak0,k1,2,,n 2.k0Ak(x)dxanb (5)求积公式的构造 第一步:找高斯点 2g(x)1,g(x)xa,g(x)xbxc,由正交性确定121)待定系数法:设0待定系数a,b,c,…..2)利用递推公式 第二步:确定求积系数Ak 1)解线性方程组 2)Ak(x)lk(x)dx,k1,2,,nab lk(x) i0iknxxi,k1,2,,nxkxi 三、本章思考题 1.插值型求积公式有何特点? 答:插值型求积公式主要用于计算定积分的值。数学推导中用拉格朗日插值函数代替被积函数,其表现形式是有限个函数值的线性组合,而组合系数恰好是拉格朗日插值基函数的定积分。(n+1)个结点的插值型求积公式的代数精度一般不超过n。用数值求积公式计算定积分可以克服牛顿—莱布尼兹公式的弱点,但是数值计算结果带有误差。在用数值求积公式设计算法时,一般要考虑到误差估计,还应该使所求的数据结果的误差得到控制。2.复化求积公式的误差是如何估计的? 答:对于复化梯形公式可根据其截断误差公式,首先求得hba,然后求nf(x)的二阶倒数,判断f(x)的二阶倒数的单调性,然后在积分区间上求得f(x)的二阶倒数的最大值就可以估计复化求积公式的误差,利用估计出的复化求积公式的误差还可以求得用复化梯形公式近似求解某一积分的有效数字有多少位。对于复化Simpson公式方法同估计复化梯形公式的误差,只是截断误差公式有所改变,此时需求出f(x)的四阶倒数然后判断其最大值。 四、本章测验题 1问题:如果用复化梯形公式计算定积分exdx,要求截断误差不超过 00.5104,试问n至少取多少? 解:复化的梯形公式的截断误差为:RTba3''hf 12RT1ba3hmaxf''(),而maxf''()max(ex)1,h 0x10x10x1n12将以上各式代入RTba3hmaxf''()可得: 0x112ba314 hmaxf''()0.51020x11212nRT解上述方程得n40.8,取n41,所以n至少取41。第二篇:数值分析学习心得体会
第三篇:清华大学数值分析实验报告
第四篇:复数值分析习题
第五篇:数值分析第六章学习小结