哈工大数值分析上机实验报告

1、哈工大数值分析上机实验报告哈工大数值分析上机实验报告52/52哈工大数值分析上机实验报告实验报告一题目：非线性方程求解纲领：非线性方程的分析解平常很难给出，所以线性方程的数值解法就特别重要。本实验采纳两种常有的求解方法二分法和Newton法及改良的Newton法。序言：（目的和意义）掌握二分法与Newton法的基本源理和应用。数学原理：对于一个非线性方程的数值解法好多。在此介绍两种最常有的方法：二分法和Newton法。对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在a,b上连续，f(a)f(b)0，且f(x)在a,b内仅有一个实根x*，取区间中点c，若，则c恰为其根，不然依据f

2、(a)f(c)5e-6);c=(a+b)/2;iff12(a)*f12(c)0;a=c;elseb=c;endR=b-a;%求出偏差k=k+1;endx=c%给出解Newton法及改良的Newton法源程序：clear输入函数f=input(请输入需要求解函数,s)%求解f(x)的导数df=diff(f);%改良常数或重根数miu=2;%初始值x0 x0=input(inputinitialvaluex0);k=0;%迭代次数max=100;%最大迭代次数R=eval(subs(f,x0,x);%求解f(x0)，以确立初值x0时否就是解while(abs(R)1e-8)x1=x0-miu*ev

3、al(subs(f,x0,x)/eval(subs(df,x0,x);R=x1-x0;x0=x1;k=k+1;if(eval(subs(f,x0,x)max;%假如迭代次数大于给定值，以为迭代不收敛，从头输入初值ss=input(mayberesultiserror,chooseanewx0,y/n?,s);ifstrcmp(ss,y)x0=input(inputinitialvaluex0);k=0;elsebreakendendendk;%给出迭代次数x=x0；%给出解结果分析和讨论：1.用二分法计算方程sinxx20在1，2内的根。(5*106,下同)2计算结果为x=；f(x)=；k=1

4、8；由f(x)知结果知足要求，但迭代次数比好多，方法收敛速度比较慢。用二分法计算方程x3x10在1，内的根。计算结果为x=；f(x)=；k=17；由f(x)知结果知足要求，但迭代次数仍是比好多。用Newton法求解以下方程a)xex100=；x计算结果为x=；f(x)=；k=4；由f(x)知结果知足要求，并且又迭代次数只有4次看出收敛速度很快。b)x3x10 x0=1；c)(x1)2(2x1)00=,x0=；x当x0=时，计算结果为x=；f(x)=；k=4；由f(x)知结果知足要求，并且又迭代次数只有4次看出收敛速度很快，实质上该方程的确有真解x=。当x0=时，计算结果为x=；f(x)=0；k

5、=9；由f(x)知结果知足要求，实质上该方程的确有真解x=，但迭代次数增添，实质被欺诈x0时，x1，就变为了方程的另一个解，这说明Newton法收敛与初值很相关系，有的时候甚至可能不收敛。4.用改良的Newton法求解，有2重根，取2(x1)2(2x1)0 x0=；并与3.中的c)比较结果。当x0=时，程序死循环，没法计算，也就是说不收敛。改1.5时，结果收敛为x=；f(x)=；k=16；明显这个结果不是很好，并且也不是收敛至方程的2重根上。当x0=时，结果收敛为x=；f(x)=；k=4；此次达到了预期的结果，这说明初值的采纳很重要，直接关系到方法的收敛性，实质上直接用Newton法，在给定相

6、同的条件和精度要求下，可得其迭代次数k=15，这说明改良后的Newton法法速度的确比较快。结论：对于二分法，只需能够保证在给定的区间内有根，使能够收敛的，当时收敛的速度和给定的区间相关，二且整体上来说速度比较慢。Newton法，收敛速度要比二分法快，可是最后其收敛的结果与初值的采纳相关，初值不相同，收敛的结果也可能不相同，也就是结果可能时时预期需要得结果。改良的Newton法求解重根问题时，假如初值不妥，可能会不收敛，这一点特别重要，自然初值适合，相同状况下其速度要比Newton法快得多。告二目：Gauss列主元消去法纲领：求解性方程的方法好多，主要分直接法和接法。本运用直接法的Guass消

7、去法,并采纳主元的方法方程行求解。序言：（目的和意）学Gauss消去法的原理。认识列主元的意。确立什么候系数要主元数学原理：因为一般性方程在使用Gauss消去法求解，从求解的程中能够看到，若a(kkk1)=0，必行行交，才能使消去程行下去。有的候即便akk(k1)0，可是其特别小，因为机器舍入差的影响，消去程也会出不定得象，致果不正确。所以有必需队列主元技，以最大可能的除去种象。一技要找行r，使得并将第r行和第k行的元素行交，以使得目前的(k1)的数比0要大的多。种列主kka元的消去法的主要步以下：消元程k=1,2,n-1,行以下步。主元，若|ark|很小，明方程的系数矩大病，出警示，提示果可

8、能不。交增广A的r，k两行的元素。arjakj(j=k,n+1)算消元aijaijaikakj/akk(i=k+1,n;j=k+1,n+1)回代程k=n,n-1,1,行以下算至此，达成了整个方程的求解。程序：本采纳Matlab的M文件写。Gauss消去法源程序：cleara=input(入系数：n)b=input(入列b：n)n=length(b);A=abx=zeros(n,1);%函数主体fork=1:n-1;%能否进行主元采纳ifabs(A(k,k)abs(t)p=r;elsep=k;endend%互换元素ifp=k;forq=k:n+1;s=A(k,q);A(k,q)=A(p,q);A

9、(p,q)=s;endendend%判断系数矩阵能否奇怪或病态特别严重ifabs(A(k,k)yipusilongdisp(矩阵奇怪，解可能不正确)end%计算消元,得三角阵forr=k+1:n;m=A(r,k)/A(k,k);forq=k:n+1;A(r,q)=A(r,q)-A(k,q)*m;endendend%求解xx(n)=A(n,n+1)/A(n,n);fork=n-1:-1:1;s=0;forr=k+1:n;s=s+A(k,r)*x(r);endt=(A(k,n+1)-s)x(k)=(A(k,n+1)-s)/A(k,k)end结果分析和讨论：26x22例：求解方程575y34。此中为

10、一小数，当105,1010,1014,1020时，321z10分别采纳列主元和不列主元的Gauss消去法求解，并比较结果。记Emax为求出的解代入方程后的最大偏差，按要求，计算结果以下：当105时，不选主元和选主元的计算结果以下，此中前一列为不选主元结果，后一列为选主元结果，下同。Emax=，0此时，因为不是很小，机器偏差就不是很大，由Emax能够看出不选主元的计算结果精度还能够够，所以此时能够考虑不选主元以减少计算量。当1010时，不选主元和选主元的计算结果以下Emax=，0此时由Emax能够看出不选主元的计算精度就不好了，偏差开始增大。当1014时，不选主元和选主元的计算结果以下00000

11、00Emax=，0此时由Emax能够看出，不选主元的结果应当能够说是不正确了，这是由机器偏差惹起的。当1020时，不选主元和选主元的计算结果以下NaN1NaN2NaN3Emax=NaN,0不选主元时，程序报错：Warning:Dividebyzero。.这是因为机器计算的最小精度为10-15，所以此时的1020就以为是0，故出现了错误现象。而选主元时则没有这种现象，并且由Emax能够看出选主元时的结果应当是精准解。结论：采纳Gauss消去法时，假如在消元时对角线上的元素素来较大（假如大于10-5），那么本方法不需要进队列主元计算，计算结果一般就能够达到要求，不然必然进队列主元这一步，以减少机器

12、偏差带来的影响，使方法得出的结果正确。告三目：Rung象生和战胜纲领：因为高次多式插不收，会生Runge象，本在出详细的例后，采用分段性插和三次条插的方法有效的战胜了一象，并且取的很好的插效果。序言：（目的和意）深刻多式插的弊端。明确插的不收性怎战胜。明确精度与点和插方法的关系。数学原理：在定n+1个点和相的函数此后结构n次的Lagrange插多式，果表明（后边的）种多式其实不是跟着次数的高升函数的迫近愈来愈好，种象就是Rung象。解决Rung象的方法平常有分段性插、三次条插等方法。分段性插：在区a,b上，定n+1个插点a=x0 x1xn=b和相的函数y0，y1，yn，求作一个插函数(x)，拥

13、有以下性：1)(xj)yj，j=0，1，n。2)(x)在每个区xi,xj上是性函数。插函数(x)称区a,b上n个数据点的分段性插函数。三次条插：定区a,b一个分划：a=x0 x1xN=b若函数S(x)足以下条件：S(x)在每个区xi,xj上是不高于3次的多式。S(x)及其2数在a,b上。称S(x)使对于分划的三次条函数。程序：本采纳Matlab的M文件写。此中待插的方程写成function的方式，以下functiony=f(x);y=1/（1+25*x*x）;写成如上形式即可，下边出主程序Lagrange插源程序：n=input(将区间分为的等份数输入：n);s=-1+2/n*0:n;%给定的

14、定点，Rf为给定的函数x=-1:1;f=0;forq=1:n+1;l=1;%求插值基函数fork=1:n+1;ifk=q;l=l.*(x-s(k)./(s(q)-s(k);elsel=l;endendf=f+Rf(s(q)*l;%求插值函数endplot(x,f,r)%作出插值函数曲线gridonholdon分段线性插值源程序clearn=input(将区间分为的等份数输入：n);s=-1+2/n*0:n;%给定的定点，Rf为给定的函数m=0;hh=;forx=-1:hh:1;ff=0;fork=1:n+1;%求插值基函数switchkcase1ifxs(n);l=(x-s(n)./(s(n+

15、1)-s(n);elsel=0;endotherwiseifx=s(k-1)&x=s(k)&xR);%精度控制j=j+1;s=0;forp=1:2(j-2);s=s+f(a+(2*p-1)*h/(2(j-1);endT(1,j)=T(1,j-1)/2+h*s/(2(j-1);%梯形公式应用form=2:j;k=(j-m+1);T(m,k)=(4(m-1)*T(m-1,k+1)-T(m-1,k)/(4(m-1)-1);endend%给出Romberg积分法的函数表I=T(m,1)结果分析和讨论：进行详细的积分时，精度取R=1e-8。1003dx。精准解I=。1.求积分x6运转程序得Romberg

16、积分法的函数表为+007*000由函数表知Romberg积分给出的结果为*107，与精准没有偏差，精度很高。2.求积分31dx。精准解I=ln3=。x运转程序得Romberg积分法的函数表为000000000000000000000从积分表中能够看出程序运转的结果为，取8位有效数字，知足要求。3.求积分1sinxdx。x直接按前面方法进行积分，会发现系统报错，出现了0为除数的现象。出现这种状况的原由就是当x=0时，被积函数分母出现了0，假如用一个适合的小数（最好不要小于程序给定的最小偏差值，可是不可以够小于机器的最大精度）来取代，能够防范这个问题。本实验取R，可得函数表为：0000000000

17、故该函数的积分为，取8位有效数字。12dx4.求积分sinx0此题的分析解很难给出，但运用Romberg积分能够很简单给出近似解，函数表为：000000000000000故该函数的积分为，取8位有效数字。结论：Romberg积分平常要求被积函数在积分区间上没有奇点。若有奇点，且奇点为第一间断点，那么采纳例3的方法，仍是能够求出来的，不然，必然采纳其他的积分方法。自然，Romberg积分的收敛速度仍是比较快的。实验报告六题目：常微分方程初值问题的数值解法纲领：本实验主要采纳经典四阶的R-K方法和四阶Adams展望-校订方法来求解常微分方程的数值解。序言：（目的和意义）经过编写程序，进行上机计算，

18、使得对常微分方程初值问题的数值解法有更深刻的理解，掌握单步法和线性多步法是怎样进行实质计算的及两类方法的合用范围和优弊端，特别是对这两类方法中最有代表性的方法:R-K方法和Adams方法及展望-校订方法有更好的理解。经过这两种方法的配合使用，掌握不相同的方法怎样配合在一同，解决实诘问题。数学原理：对于一阶常微分方程初值问题dy/dxf(x,y)y(x0)（1）y0的数值解法是近似计算中很中的一部分。常微分方程的数值解法平常就是给出定义域上的n个等距节点，求出所对应的函数值yn。平常其数值解法可分为两大类：单步法：这种方法在计算yn+1的值时，只需要知道xn+1、xn和yn即可，即可算出。这种方

19、法典型有欧拉法和R-K法。多步法：这种方法在计算yn+1的值时，除了需要知道xn+1、xn和yn值外，还需要知道前k步的值。典型的方法如Adams法。经典的R-K法是一个四阶的方法。它最大的长处就是它是单步法，精度高，计算过程便于改变步长。其弊端也很明显，计算量大，每行进一步就要计算四次函数值f。它的详细的计算公式如式（2）所示。四阶Adams展望-校订方法是一个线性多步法，它是由Adams显式公式yn1yn(K12K22K3K4)h/6K1f(xn,yn)K2f(xnh/2,ynK1h/2)（2）K3f(xnh/2,ynK2h/2)K4f(xnh,ynK3h)和隐式公式构成，其计算公式如式（

20、3）所示。展望yn1yn(55fn59fn137fn29fn3)h/24（3a）求导fn1f(xn1,yn1)（3b）校订yn1yn(9fn119fn5fn1fn2)h/24（3c）求导fn1f(xn1,yn1)（3d）将局部截断偏差用展望值和校订当来表示，在展望和校订的公式中分别以它们各自的阶段偏差来进行填补，可希望的到精度更高的修正的展望-校订公式为：展望pn1yn(55fn59fn137fn29fn3)h/24修正mn1pn1(cnpn)251/270求导fn1f(xn1,mn1)校订cn1yn(9fn119fn5fn1fn2)h/24修正yn1cn1(cn1pn1)19/270求导fn

21、1f(xn1,yn1)因为开始时无展望值和校订当能够利用，故令p0=c0=0，此后按上边进行计算。此方法的长处是能够减少计算量；弊端是它不是自开始的，需要先知道前面的四个点的值y0,y1,y2,y3，所以不可以够独自使用。其他，它也不便于改变步长。程序设计：本实验采纳Matlab的M文件编写。此中待求的微分方程写成function的方式，以下functionyy=g(x,y);yy=-x*x-y*y;写成如上形式即可，下边给出主程序。经典四阶的R-K方法源程序clear%步长采纳h=;%初始条件，即x=0时，y=1。y(1)=1;%求解区间a=0;b=2;%迭代公式forx=a:h:b-h;k

22、1=g(x,y(x-a)/h+1);%g(x)y(x)，下同。k2=g(x+h/2,y(x-a)/h+1)+h/2*k1);k3=g(x+h/2,y(x-a)/h+1)+h/2*k2);k4=g(x+h,y(x-a)/h+1)+h*k3);y(x-a)/h+2)=y(x-a)/h+1)+h*(k1+2*k2+2*k3+k4)/6;end四阶Adams展望-校订方法源程序%步长采纳h=;%初始条件y(1)=1;%求解区间a=0;b=2;%应用RK迭代公式计算初始值y0,y1,y2,y3forx=a:h:a+2*h;k1=g(x,y(x-a)/h+1);k2=g(x+h/2,y(x-a)/h+1)+h/2*k1);k3=g(x+h/2,y(x-a)/h+1)+h/2*k2);k4=g(x+h,y(x-a)/h+1)+h*k3);y(x-a)/h+2)=y(x-a)/h+1)+h*(k1+2*k2+2*k3+k4)/6;end%应用展望校订法求解c(4)=0;%校订初值p(4)=0;%展望初值f(1)=g(a+0*h,y(1);g(x)y(x)，且将该值存在数组f中。f(2)=g(a+1*h,y(2);f