数值分析各算法流程图

数值分析各算法图 数值分析各算法流程图 一、插值 1、 拉格朗日插值流程图:( 相应程序:lagrintp(x,y,xx)) 输入及xyinx,(1,,),？ii y,0 kn,1,,？ P,1 jn,1,,？ 是 kj,? 否 PPxxxx,,,*()/()jkj yyPy,，*k 输出xy, 2、 牛顿插值流程图 (1)产生差商的算法流程图(相应程序:divdiff(x,y)) 输入xyin,,(1,,),？ii Diyin(,1),1,2,,,,？i jn,2,,？ ijn,,,？ DijDij(,1)(1,1),,,, Dij(,),xx,iij,，1 输出矩阵D 注:1、另一程序divdiff1(x,y)，输出的矩阵包含了节点向量。而divdiff(x,y)不含节点向量。 2、另一程序tableofdd(x,y,m)，输出的是表格形式，添加了表头。 (2)非等距节点的牛顿插值流程图(相应程序:newtint11(x,y,xx,m)) 、 输入及xyinmx,,(1,,),,,？ii 调用差商表程序，产生矩阵D sDw,,(1,1),1 im,1,,？ wwxx,,*(), i ssDiiw,，，，(1,1)*, 输出Ds, 注:1、虽然程序newtint11(x,y,xx,m)考虑了多种情形，看上去很复杂，但基本流程结构还是如上图所示。 2、程序中调用的子程序是divdiff。 若调用的子程序是divdiff1的话，流程图中的第三，第四，第五步要相应的改一下数字。 (3)求差分表的流程图(相应程序:difference(y,m)) 输入及yinm,(1,2,,),？i Diyin(,1),1,2,,,,？i jm,，2,3,,1？ innj,,,1,,？ DijDijDij(,)(,1)(1,1),,,,, 输出矩阵D 注:1、difference输出的是矩阵D。而另一程序tableofd(y,m)，输出的是带有表头的差分表。 (4)牛顿向前插值的流程图(相应程序:newtforward1(x,y,xx,m)) 输入xyinxm,,(1,,),,,？ii xx,1 hxxt,,,,21h 是 tt,,01?或输出出错信息 否 执行求差分的子程序，并将输出结果赋值给矩阵D sDw,,(1,1)1， im,1,,？ wwtii,,，*(1)/ sswDii,，，，*(1,1) 输出s 注:1、程序newtforward1(x,y,xx,m))的结构与上述流程图一致，xx可以是数组。 2、另一程序newtforward(x,y,xx,m))先求出插值多项式，再求插值多项式在插值点的函数值。基本结构还是和上面的流程图一样。 (5)牛顿向后插值的流程图(相应程序:newtbackward1(x,y,xx,m)) 输入xyinxm,,(1,,),,,？ii xx, 1hxxt,,,,21h 是 tt,,01,或?输出失败信息 否 执行求差分的子程序，并将输出结果赋值给矩阵D sDnw,,(,1)1， im,1,,？ wwtii,，,*(1)/ sswDni,，，*(,1) 输出s 注:1、程序newtbackward1(x,y,xx,m))的结构与上述流程图一致，xx可以是数组。 2、另一程序newtbackward(x,y,xx,m))先求出插值多项式，再求插值多项式在插值点的函数值。基本结构还是和上面的流程图一样。 3、Hermite插值流程图 (1) 已知条件中一阶导数的个数与插值节点的个数相等时的Hermite插值流程图。(相应程序: hermiteint(x,y,y1,xx)) ,输入及xyyinx,,,(1,2,,),？iii y,0 in,1,2,,？ lm,,1,0 kn,1,2,,？ 否 ki,? 是 llxxxx,,,*()/()kik mmxx,，,1/()ik 2,yyxxmyxxyl,，，,,，，,12*****，，，， ,,iiii，， 输出y 注:hermiteint(x,y,y1,xx)程序中的基本结构与上述流程图一样，但还考虑了不同情形下的输出 结果:若输入前3个参数，则输出的是插值多项式，若输入了4个参数，则输出xx相应的插 值函数值。 (2) 已知条件中一阶导数的个数小于或等于插值节点的个数时的Hermite插值流程图。(相应程 序:HermiteInt1(x,y,x1,y1,xx)) 输入及xyinxykmx,,(1,2,,),1,1,(1,2,,),,？？iikk y,0 in,1,2,,？ I,0 km,1,2,？ 否 xx,1?ik 是 Ik, 是否I,0? nm xx,1xx,jkyyy,，,,,,i 1xxxx,,,,11jkijik,ji ，，,,nmnmxxxx,, 11ij,,k,,1yyxxy，，,，,,,，,,,,,,ii,,,, xxxxxxxx,,,,,,11,,11jkjkijiiijii,,kk,,,,,,jiiijiiikk,,，， km,1,2,,？ nm xx,xx,ijl1yyxxy,，,,,,，，,,ikk xxxx,,,,11jlijiikkl,,jilkk 输出 y 注:上面的是基本流程图，事实上，HermiteInt1的详细流程图如下所示 输入xyin,,(1,2,,),,？ii xykmx1,1,(1,2,,),？及kk y,0否I,0? 是 in,1,2,,？t,0 I,0 yystxxy,，,,,,,(1()) ii km,1,2,？ 否 xx,1?ikkm,1,2,,？ 是 Ik,r,1 in,1,2,,？ 是st,,1,0xx,1?ik 否 in,1,2,,？rrxxxx,,,,()/(1)iki 是ji,? 否 jm,1,2,,？ssxxxx,,,,()/()jij 是 xx11?,ttxx,，,1/()jkij 否 rrxxxx,,,,(1)/(11)jkj km,1,2,,？ 是 yyrxxy,，,,,(1)1kI,?kk 否 ssxxxx,,,,(1)/(1)kik 输出 yttxx,，,1/(1)ik 4、分段低次插值流程图 (1)分段线性插值流程图(相应程序:fendlineint(x,y,xx)) 输入及xyinx,(1,,),？ii iaibn,,1, 否 ibia,,1? 是 iaib，im,的整数部分 执行求拉格郞日插值的 2 ,程序，计算以为节iaib 点的线性插值函数在插否ibim,xx,?im ()值点的函数值Px1 是 iaim,输出Px()1 注:1、fendlineint(x,y,xx)的流程结构与流程图一致，参数xx表示插值点，只能是一个数。 2、另一程序PiecewLinInt(x,y,xx)，可以输出每个子区间上的插值函数表达式。 这里的xx可以是一个数组。 此程序调用了程序fendlineint(x,y,xx)。 3、程序PicecLPlot(f,x)是画出原函数与插值函数在插值区间上的图像。 (2)分段三次Hermite插值的流程图(相应程序:PieceThHermInt(x,y,y1,xx)) 输入及xyinx,(1,,),？ii iaibn,,1, 否ibia,,1? 是 iaib， 执行求插值的程序，Hermiteim,的整数部分2 计算以为节点的次iaibHermite,3 插值函数在插值点的函数值xHx()3否 ibim,xx,?im 是 输出Hx()3iaim, 注:1、程序PieceThHermInt(x,y,y1,xx)调用了程序hermiteint(x,y,y1,xx)， 2、参数xx可以是数组。 5、三次样条插值流程图 (1)第一种边界条件下的三次样条插值流程图(相应程序:cubicspline12) ,,输入及xyinyyx,(1,),,,？iin1 hxx,,iii，1 1,,1in,,？[,]()/fxxyyh,,iiiii，，11 ,gfxxyh,,6[,]/，，11211 ,gyfxxh,,6[,]/，，nnnnn,,11 h i,,,1i，hh ，1ii ,,,,12,,1？in,,,,11ii 6,,[,][,]gfxxfxx，，，,11iiiii ，hh，1ii g21，，，，1 ,,,,g,,2211 ,,,, ,,,,Dg,,,？？？？ ,,,,g2,, n,1nn,,22,,,, ,,,,g12 n，，，， ,1MDg, kn,,2,,1？ 是Ik,xx,?k 否 In, 33()()xxxxMxxMxx,,,,,,,,22IIIIII,,,111yMMyhyh,,，,，,,，,, IIIIII,,,,1111,,,,6666hhhh,,,,IIII,,,,1111 输出xy, (2)第二种边界条件下的三次样条插值流程图(相应程序:cubicspline22) ,,,,输入及xyinyynx,(1,),,,,？iin1 hxx,,iii，1 1,,1in,,？[,]()/fxxyyh,, iiiii，，11 hi,,i ，hh，1ii 11,,2,,,,？in,,ii 6[,][,],,gfxxfxx，， ，，，121iiiii，hh，1ii ,,2,,gy,，，，，1111,,,, 2g,,222,,,, ,,,,Dg,,？？？？,,,,, 2g,,nnn,,,333,,,, ,,,,,,2gy,nnnn,,,222，，，，,, ,1 MDg, kn,,2,,1？ 是 xx,?Ik,k 否 In, 33 ()()xxxxMxxMxx,,,,,,,,22IIIIII,,,111yMMyhyh,,，,，,,，,,IIIIII,,,,1111,,,, 6666hhhh,,,,IIII,,,,1111 输出 y 3)第三种边界条件下的三次样条插值流程图(相应程序:cubicspline32) ,,,,输入及xyinyyx,(1,),,,？iin1 hxx,,iii，11,,1in,,？ [,]()/fxxyyh,,iiiii，，11 hi ,,i，hh，1ii 11,,2,,,,？in,,ii 6 [,][,],,gfxxfxx，，，，，121iiiii，hh，1ii hn,1,,,,,,,1nnn,,,111 hh，11n, 6gfxxfxx,,[,][,]，，nnn,,1121 hh，11n, 2,g，，，，11 ,,,,2g,,222,,,, ,,,,Dg,,？？？？, ,,,,2g,,nnn,,,222,,,, ,,,,2gnn,,11，，，，, ,1 MDg, kn,,2,,1？ 是Ik, xx,?k 否 In, 33 ()()xxxxMxxMxx,,,,,,,,22IIIIII,,,111yMMyhyh,,，,，,,，,,IIIIII,,,,1111,,,, 6666hhhh,,,,IIII,,,,1111 输出 y (4)三次样条插值总程序流程图(相应程序:cubicspline) 输入数据flag 是 flag,1?执行子程序“第一类边界条件下的样条三次插值” 否 是flag,2?执行子程序“第二类边界条件下的样条三次插值” 否 是flag,3?执行子程序“第三类边界条件下的样条三次插值” 否 输出失败信息 二、拟合 1、用最小二乘法求m次多项式拟合的流程图(相应程序:lsappro1) 输入及 ,(1,2,,),xyinmx,？ii n，,2ij 求法矩阵:Aijximjm(,),1,,1;1,,,1,,，,，？？,k,1k n ,1i求常数项:bixyim(),1,,1,,，？,kk ,1k ,1求回归系数:aAb, m，1 k,1yaxinaaa1,1,2,,.,,,,,？？其中，，,ikim11， k,1 n 2erroryy,,(1),ii ,1i 输出 ,aerror 注:另一程序lsappro较复杂，输出的是拟合多项式的表达式，误差，和已知数据及拟合多项 式的图像。 2、用最小二乘法求线性模型拟合的流程图(相应程序:lsmethod1) 输入fimxyjn(1,),,(1,),,？？ijj aijfximjn(,)(),1,;1,,,,,？？ij ,,Aaabay,,,,, ,1 cAb,, m yjfxjn1()(),1,,,？ ,iji,1 n2erroryy,,(1) ,jj,j1 输出 ,cerror 注:lsmethod2程序中还输出了拟合多项式的表达式，图形，及在插值点的函数值。 abx，(3)求最小二乘法求形如的拟合曲线(相应程序:leexp(x,y)) ye, 输入xyjn,(1,),？jj 执行用最小二乘法求次1 拟合多项式程序，将法方 程组的解赋值给向量C CCx(1)(2)，,iyejn1,1,,,？ j n2erroryy,,(1) ,jj,j1 输出 ,Cerror 注:leexp(x,y)程序输出的是拟合函数的表达式，还可以自由的输出拟合函数与已知数据的图 形，以及可根据需要输入插值点的横坐标，输出相应的拟合函数值。 注2: blappro是最佳一致逼近的matlab程序，流程图没写。 三、数值积分 1、复合Newton-Cotes求积的流程图(相应程序:NewtonCotes1) 输入 ,,,,fabnm hban,,()/ xaihinn,，,,，(1),1,,,1？i in,1,,？ xx,ii，11(1),xxj,，,, jim1,,1jm,，？ 1(1),yfx,jj m，1 yhyc,,1,ikmkk,1 n Iy,,i i,1 输出I 2、变步长梯形求积公式(相应程序:ComTrap(fun,a,b,delta,m)) ComTrap1 输入fabm,,,,, hba1,, xaxb1,1,,12 yfxi1(),1,2,,ii 2h1,,ty1 ,i12,i1 error,1 im,1,,？ 否 输出Ierror,3?, 是 hh21/2, xxh212,，jj,1ij,1,,2？ yfx2(2),jj i,12 tthy,，,/222,iij，1 j,1 errortthh,,,,12ii，1 ixxk11,1,,2,,？21kk, i,1 xxk12,1,,2,,？2kk 输出失败信息 3、龙贝格算法流程图(相应程序:romberg(fun,a,b,delta,m)) 输入fabm,,,,, hba1,, xaxb1,1,,12 yfxi1(),1,2,,ii 2 h1,,Iy1,i11 2,i1 error,1 im,1,,？ 否输出Ierror,,? 是 hh21/2, xxh212,，jj ,1ij,1,,2？ yfx2(2),jj i,12 IIhy,，,/222,iij，1,1,1j,1 ji,，2,,1？errorIIhh,,,,12iiii，，1,1, j,1 4,,IIijij，,,1,1,1iI,xxk11,1,,2,,？ij，1,j,121kk, 41, i,1 xxk12,1,,2,,？2kk 输出失败信息 四、非线性方程的数值解法 1、二分法的流程图(相应程序:bisection) 输入fabn,,,,,,,12 kn,1,,？ midab,，()/2 errorab,,()/2 是 errorfmid,,,,或()?输出mid12 否 是fmid()0?,输出mid 否 是 fafmid()()0,,bmid, 否 amid, 输出失败信息 注:1、程序bisection求确定函数fun在区间[a,b]内的一个近似根 2、程序bisection1(fun,a,b,delta1,delta2,n)是将逐步搜索法与二分法结合起来，确定函数fun 在区间[a,b]内的多个近似根，此程序调用了两个子程序ssf和bisection。 2、定点迭代的流程图(相应程序:fpiteration) 输入fxn,,,,0 kn,1,,？ xfx,()0 errorxx,,0 是error,,?输出x 否 xx,0 输出失败信息 3、 牛顿迭代的流程图(相应程序:newtiteration) , 输入ffxn,,,,,0 kn,1,,？ 是 ,fx()0?,输出出错信息0 否 ,xxfxfx,,()/() 000 errorxx,,0 是 error,,?输出x 否 xx,0 输出失败信息 4、Aitken加速法的流程图(相应程序:AitkenM) 输入fxn,,,,0 in,1,,？ yfxzfy,,(),()0 2 ()yx,0xx,,i0 zyx,，20 errorxx,,i0 是 输出xerror,,?i 否 xx,0i 输出失败信息 5、牛顿下法的流程图(相应程序:NewDecentiteration) ,输入ffxn,,,,,,,0 kn,1,2,,？ t,1 是,fx()0?,输出出错信息0 否 fx()0 xxt,,,0,fx()0 否否t,,?fxfx()()?,tt,/2 0 是 输出出错信息 是 是 输出xxx,,,,0 否 xx,0 输出失败信息 6、正割法的流程图(相应程序:SecantIteration) 输入fabn,,,,,,,12 kn,2,,？ ()ba, xafa,,()*fbfa()(), 是输出xxbfx,,,,,或()?12 否 是 fafx()()0,,bx, 否 ax, 输出失败信息 五、线性方程组的数值解法 1、Gauss顺序消元法的流程图(相应程序:gauss) 输入ab,，，，，ijin nn， kn,,1,,1？ 是,10 a,10?输出失败信息kk 否 aa,kjik,1,,,aaikn,,,，？ijij akk 1,,jkn,，？ ba,kik,1,,bbikn,,,，？ii akk xba,/ nnnn knn,,,1,2,,1？ n bax,,kkii ik,，1x,ka kk 输出xkn,1,,,？k 2、利用Doolittle分解求线性方程组解的流程图 (1)矩阵Doolittle分解的流程图(相应程序:DoolFactorization) 输入a，，ij，nn cajn,,,1,？11jj BI, nn， ai1,2,,bin,,？ i1c11 kn,2,,？ k,1 cabcjkn,,,,,,？ ,kjkjkiij,,i,1 k,1 abc,,ikijjk j,1,1,,？bikn,,，ik ckk 输出　bc, ，，，，ijij，，nnnn (2)求解系数矩阵为下三角矩阵的线性方程组的流程图(相应程序:lowerTriangular) 输入ab,，，，，ijin，nn 输出失败信息:是a,0? 11系数矩阵不可逆 否 xba,/1111 kn,2,,？ 是输出失败信息:a,0?kk 系数矩阵不可逆 否 k,1 bax,,kkjj j,1x,k akk 输出xkn,1,,,？k (3)求解系数矩阵为上三角矩阵的线性方程组的流程图(相应程序:upperTriangular) 输入ab,，， ，，ijin，nn 输出失败信息:是a,0?nn 系数矩阵不可逆 否 xba,/nnnn knn,,,1,2,,1？ 是输出失败信息: a,0?kk 系数矩阵不可逆 否 n bax,,kkjjjk,，1 x,ka kk 输出xkn,1,,,？k (4)利用Doolittle分解求线性方程组解的主程序流程图(相应程序:TriangularDec) 输入ab,，，，，ijin，nn 执行分解子程序，并将下三角矩Doolittle 阵赋值给矩阵将上三角矩阵赋值给矩阵BC, 执行求角系数矩阵为下三角的线性方程组 的子程序，求,的解，并将解赋值给BXby 执行求系数矩阵为上三角的线性方程组解 的子程序，求,的解，并将解赋值给CXyx 输出x 3、 追赶法的流程图(相应程序:zgTridiagonal) 输入abcd,,,，，，，，，，，iiii,,11nnnn kn,2,,？ 是b,0?输出失败信息k,1 否 cakk,,11,,,bbkk bk,1 dakk,,11 ,,ddkkbk,1 是 b,0?输出失败信息n 否 xdb,/ nnn knn,,,1,2,,1？ dcx, kkk，1x,kb k 输出xkn,1,,,？k 4、列主元消去法的流程图(相应程序:MColumnPivoting1) 输入ab,，，，，ijinnn， kn,,1,,1？ dk, ikn,，1,,？ 否aa,? ikdk 是 w di, caeb,,jkjk ,1,,,aajnbb,,,？ kjdjkd acbe,,djjd 是,10 a,10?输出失败信息kk 否 aa, kjikxba,/,1,,,aaikn,,,，？nnnnijij akk 1,,jkn,，？knn,,,1,2,,1？ ba, kikn,1,,bbikn,,,，？iibax, a,kkiikkik,，1x, kakk 输出xkn,1,,,？k 5、平方根法的流程图 (1)Cholesky分解的流程图(相应程序:CholeskyFactorization) 输入　a，，ij，nn cinjn,,,0,1,,;1,,　？？ij in,1,,？ jin,,,？ 是 j,1?sa,11 否 i,1 sacc,,,ijkikjk,1 否 csc,/ji,?ijii 是 是s,0?输出失败信息 否 cs,ii 输出　c，，ijnn， (2)平方根法主程序的流程图(相应程序:CholeskyMethod) 输入ab,，，，，ijin，nn 执行分解子程序，并将Cholesky 输出的上三角矩阵赋值给矩阵B 执行求角系数矩阵为下三角的线性方程组 ,的子程序，求,的解，并将解赋值给BXby 执行求系数矩阵为上三角的线性方程组解 的子程序，求,的解，并将解赋值给BXyx 输出x 六、解线性方程组的迭代法 1、Jacobi迭代的流程图(相应程序:JacobiIteration) 输入(),(),0,,maxabx,ijnnin， kn,1,,？ 是 ,10a,10?输出失败信息kk 否 k,1,,max？ n xbaxain,,,(0)/,1,,？,iiijjii ,1j,ji 是max0?xx,,,输出x，，iiin 1,,in 否 xxin0,1,,,,？ ii 输出失败信息 注:另一程序JacobiIteration1是根据Jacobi迭代的定义,按矩阵运算方式来编写。 2、Seidel迭代的流程图(相应程序:GaussSeidel) 输入(),(),0,,maxabx,ijnnin， kn,1,,？ 是 ,10a,10?输出失败信息kk 否 k,1,,max？ axin,,0,1,,？ ii n xbaxain0(0)/,1,,,,,？ ,iiijjii,1j,ji 是 max0?ax,,,输出x0，，iiin1,,in 否 输出失败信息 注: 另一程序GaussSeidel 1是根据Gauss Seidel迭代的定义,按矩阵运算方式来编写。 3、超松弛迭代的流程图(相应程序:SorM) 输入(),(),0,,,maxabx,, ijnnin， 0a0,,aa？,,,,,,11121n ,,,,,,,,,,,a0,,a ？？？2122,,,,,,DLU,,,,, ,,,,,,？？？？0,a,nn1,,,,,,, ,,,,,,,,aa？0a0,nnn1,1nn,,,,,, ,1BDLUD,,,，,()[(1)],,, ,1fDLb,,,() , k,1,,max？ xBxf,,，0 是Axb,,,?输出x , 否 xx0, 输出失败信息 注:此流程图比较简单，在程序中SorM中还考虑了一些可能出错的情形，如参数,应该在(0,2) 之间等。 七、常微分方程求解 1、欧拉法的流程图(相应程序:euler) 输入fxybn,0,0,, hbxn,,(0)/ xxyy,,0,011 kn,1,？ xxh,，kk，1 (,)yyhfxy,，,kkkk，1 输出xykn,,1,1.,，？ kk 2、改进欧拉法的流程图(相应程序:EulerModi) 输入fxybn,0,0,, hbxn,,(0)/ xxyy,,0,011 kn,1,？ xxh,，kk，1 tyhfxy,，, (,)kkk hyyfxyfxt,，，(,)(,)，， kkkkk，，112 输出xykn,,1,1.,，？ kk 3、四阶的龙格,库塔法的流程图(相应程序:RungeKutta4) 输入fxybnc,0,0,,,() i13 hbxn,,(0)/ xxyy,,0,011 kn,1,？ Kfxy1(,),kk KfxchychK2(,1),，,，,,56kk KfxchychKchK3(,(12)),，,，,,，,,789kk KfxchychKchKchK3(,(123)),，,，,,，,,，,, 10111213kk yyhcKcKcKcK,，,，,，,，,*1234，，，11234kk xxh,， ，1kk 输出xykn,,1,1.,，？kk 注:1、程序RungeKutta2是2阶龙格库塔方法求解一阶常微分方程。 2、两个程序都有画图部分，但没有计算误差。