2016春计算方法实验指导-作业

/计算方法实验指导材料目录TOC\ｏ"1－3"\ｈ\z\ｕHYPEＲＬＩNK\l＂_Toc４5161６0１０”第1次解线性方程组的直接解法ﻩPAＧERＥＦ＿Toc4５16１6０1０\h—2-HYPEＲLINK1．1例题 PＡGEREＦ_Toc451616011\h—2-ＨYPＥRLIＮK\l”_Ｔoｃ451616０12"1.2Matlab解线性方程组常用命令介绍 ＰＡGEREF_Tｏc451616０１２\h－6-HYPEＲＬINK＼l"_Toc４51６16０13”1。3实验作业 ＰAＧEＲEF_Toc4５16１6013\h—７－HYＰＥRＬＩNK\l”_Toc45１6１6０１４＂第2次解线性方程组的迭代解法ﻩPAGERＥＦ_Tｏc4５16１6014\h－8—HYPERＬINＫ\l＂_Toｃ451616015"2。1例题ﻩPAＧＥＲEＦ＿Toc451６1６０15\h－８—HＹPERＬINK＼l"_Toｃ４５１616016"2.2Matlaｂ迭代解法常用函数介绍 ＰＡGＥＲＥF＿Ｔｏｃ4５１6１6０16\h-12—HYPERLINK\l＂_Tｏc4５161601７"2．３实验作业ﻩＰAGＥREF_Toc４51616０1７\h-１2－HYPERLINK\ｌ”_Toｃ4５１61６018”第3次非线性方程求根 PＡGＥRＥF_Toc4５1616018＼h—13-HYPＥRLIＮK＼ｌ"_Ｔoｃ451616０１9”３.１例题ﻩPAGEＲＥF_Tｏc45１６1６0１９＼h－1３-ＨYＰERLINK\ｌ”_Toc4５１6１6020”3。2Mａｔlab非线性方程求根的命令ﻩPＡGEREF＿Toｃ45161６02０\h-22－HYPERLINK\ｌ＂_Ｔoc45161６０21＂3。3实验作业 PAGEREＦ_Toc45１6１602１\ｈ—22-ＨYPＥRＬＩNK＼ｌ＂_Toｃ４516１602２＂第4次插值法ﻩPAGＥREF＿Toc4５161６02２\h２3HYPＥRLINK\l＂＿Ｔｏc4516160２3”4。１例题ﻩＰＡGEREF_Toｃ451616０２3＼h2３HYPERLIＮＫ\ｌ"_Toｃ4５1６16024"4。２Ｍａｔlaｂ插值函数介绍 ＰAGEＲＥＦ_Tｏc4516１6０２4＼h２9HYＰEＲLIＮK\l”＿Toｃ4５161602５"4。3实验作业ﻩＰAGEREF_Toc45１6１6025\h30ＨYPＥRLINＫ＼l＂＿Ｔoc4516160２６＂第5次利用最小2乘法进行曲线拟合ﻩPAGEＲEＦ＿Toc4５1６160２6\h３1HYPERＬINＫ\l"_Toｃ451616027”5.1例题ﻩPAGEREF_Tｏc45161６０27＼ｈ31HYＰERＬINＫ＼l"＿Tｏc45161６0２8＂5。2Matｌab数据拟合命令介绍ﻩPAGＥRＥＦ＿Toc451６16０28\h34ＨYＰＥＲLINK\l”_Tｏｃ４５16１60２９＂5．３实验作业 ＰＡGEREF_Tｏｃ451616029＼h35HＹPEＲLINK\l＂＿Toｃ４51616０30＂第6次数值积分与数值微分 PAＧEＲEF_Toc4５1６１6030\ｈ36HYＰERLINＫ\l"＿Toc45１616０３1”６．1例题ﻩPＡGEREF_Tｏｃ451６１６0３1\h３6HＹPERLIＮＫ\l”_Toc4５1６1６0３2”６。2Mａtlab数值积分函数介绍 PＡGEREF_Toｃ451616032＼h42HＹＰERLINＫ\l"_Toc４51616033＂6.３实验作业ﻩPAGEREＦ_Tｏｃ45１６16033\ｈ4３第１次解线性方程组的直接解法Gａuss消元法算法如下:Gausｓ列主元消元法算法:在以上的算法中，增加调换行的步骤即可。在每步消元之前要按列选出主元。1。1例题例1.1用Gauss消元法解方程组:解：直接建立求解该方程组的M文件Gausｓ．m如下％求解例题１．1%高斯法求解线性方程组Ax=b％Ａ为输入矩阵系数，b为方程组右端系数％方程组的解保存在x变量中％先输入方程系数Ａ=［12３；275；1４９];b=［16—３］’；[m，ｎ］=siｚe（A）；%检查系数正确性iｆm~=nerroｒ('矩阵Ａ的行数和列数必须相同'）；rｅturｎ；enｄifm～＝ｓize(ｂ)ｅｒror（'b的大小必须和A的行数或A的列数相同’);return;end％再检查方程是否存在唯一解ifrａnk(A)~=rａnｋ([A,b])eｒｒｏr（＇A矩阵的秩和增广矩阵的秩不相同，方程不存在唯一解＇)；rｅｔurｎ;eｎｄ%这里采用增广矩阵行变换的方式求解ｃ＝n+1；A（：,ｃ）＝b;%％消元过程fork=1：ｎ－1A（ｋ+1：n，ｋ:c)=Ａ(k+1:n，k:c）－（A（ｋ+1:n，k)/A(k,k））*A（k，k：c）；Eｎｄ％%回代结果x=zerｏs（length（b)，1）;x(n）=A（n，c）/A(ｎ,n）;fｏrk＝n－1:－1：1x（k)＝(A(k，c）-Ａ（ｋ,ｋ+1：n)*x(k+1：n））/A(k，ｋ);eｎd％显示计算结果ｄiｓｐ(＇ｘ=＇）;diｓp（x）；直接运行上面的M文件或在Matｌab命令窗口中直接输入Gａuss即可得出结果。在Ｍaｔlａb命令窗口中输入Gauss得出结果如下:〉〉Gauｓsx=2.0０001.0000－1。0０００扩展:Matｌａb求解线性方程的几种命令如下(方程组的一般形式可用矩阵和向量表示成，但运用下列方法的前提必须保证所求解的方程为恰定方程，即方程组存在唯一解)。运用求逆思想：：或;左除法：原理上是运用高斯消元法求解，但Matlａb在实际执行过程中是通过分解法进行的(即先将矩阵A作分解,再回代计算）：；符号矩阵法:这种计算方法最接近精确值，但计算速度最慢：;将矩阵施行初等行变换化成行简化阶梯形的办法：可以这样实现之:；.上面四种常用的办法示例如下：＞>Ａ=[１２3；２75；149］%上面示例方程组系数Ａ=1２３275１4９〉〉b=[16－3]'%方程组右端的系数b=1６－3>>x1_1＝iｎv(Ａ）*ｂ，x1＿2＝A^(—１）＊ｂ%方法一，求逆思想ｘ1＿1=2.0００01.０00０—1.０000x１_2=2．００001。0０0０-1.０0０0＞>ｘ2＝A\b％方法二，左除思想x2＝21-1〉〉x3=sｙｍ（A）\sym（b)%方法三，符号法ｘ３=２1—1>〉C=[Ａ，b］，rｒｅｆ(C）％方法四,行简化阶梯形思想，最后输出结果的一列为解Ｃ=1２3127561４9－3ans=10０2０1０1001－１例1。2，求．解:输入矩阵：>＞ｘ＝[１0。５－０。3］ｘ＝1.００000。5000-0。３000计算其1-范数〉＞norm_1＝norｍ(x,1）ｎoｒｍ_１=１。8０00计算其2—范数>〉noｒm_2＝norm（x）nｏrm_2=１.1576计算其无穷大范数>>noｒm_ｉnf=noｒm(x,iｎf）norm_iｎf=1还可以计算其无穷小范数(即各分量绝对值中的最小值)>〉nｏrm_ｍinuｓInf=norm(ｘ,-inf）ｎorm_minusInf=0．3０00例1。3设矩阵，求。解:先输入矩阵:〉〉A＝［-210；1-21；01—２］A＝-2101-2１01－２求Ａ的１－范数(列和范数):〉＞noｒｍ_1=nｏrm（A,1）norｍ＿1=4求解A的无穷大范数(行和范数)：>>nｏrm_iｎf=noｒｍ(A，inｆ)nｏrm_iｎf=4求A的2－范数（最大特征值）:>〉norm_２=ｎｏｒm（A,2）norm＿２=3。4142还可以求解A的Ｆ—范数：>〉ｎｏrm_Ｆ=norm（A,'fro’）norm_F=4谱半径可以通过按其概念进行计算：对其特征值的绝对值取最大值即可：〉＞R_A=max(ａbs(eig（Ａ)））R_A＝3.4１421。2Matlaｂ解线性方程组常用命令介绍１．求秩命令rank在解线性方程组时，为了判断是否存在解，我们应先判断矩阵的秩。调用格式为:c=raｎk（A)2.矩阵零空间向量null当线性方程组的秩时，方程组有无穷多个解,使用x=null（A）可得到满足的一个解向量，可为符号矩阵或数值矩阵:x=nｕll(A)或x＝ｎｕll(syｍ（A))3.方阵的LU分解命令lu调用格式为：［L，U］=lu（A）L为准下三角矩阵，U为上三角矩阵，满足A=LU.［L，Ｕ,Ｐ］＝ｌu（A）L为下三角矩阵，U为上三角矩阵，Ｐ为变换方阵元素位置的换位阵,满足PA=PＬ．其它调用格式请用ｈｅlpｌｕ获得更多信息.４.Chｏlsky分解chｏlL=chol（Ａ)其中L为一个下三角矩阵，满足.必须为正定矩阵．5。条件数condc＝cond（A，p)A可为向量或矩阵，P取值为下列之一:１——向量或矩阵的返回１—条件数.2——返回向量或矩阵的2—条件数。inf——返回向量或矩阵的—条件数．’fro'--返回向量或矩阵的Ｆ—条件数。6。奇异值分解svｄ［U,S,Ｖ］＝svd（A）将A分解为正交矩阵U,对角矩阵Ｓ和正交矩阵Ｖ的乘积，使得A=USVＴ.7．线性方程组的符号解ｌiｎsｏｌveX＝linsolve(A，b）它等价于X=ｓym(Ａ)＼sｙｍ（b）。返回方程组的符号解。1。3实验作业实验名称:解线性方程组的列主元素高斯消去法．实验目的:通过数值实验，从中体会解线性方程组选主元的必要性，以及方程组系数矩阵和右端向量的微小变化对解向量的影响.实验内容：解线性方程组：。实验要求:(1）用Mａｔlａｂ语言编写程序，用列主元高斯消去法求解上述方程组，输出Ａx=b中矩阵A、b、解向量ｘ及ｄetA．（2）将方程组中系数3.01改为３.０0,０。987改为０．990，用列主元高斯消去法求解系数变换后的方程组，输出列主元行交换次序、解向量x及dｅｔA,并与=1\*GB2（1）中结果比较.(３)用MAＴLAB的内部函数inｖ求出系数矩阵的逆矩阵，再输入命令ｘ=inv(A)*b,即可求出上述各个方程组的解，并与列主元高斯消去法求出的解进行比较，体会选主元的方法具有良好的数值稳定性．用MATLAＢ的内部函数dｅt求出系数行列式的值,并与(１）、(2)、中输出的系数行列式的值进行比较．

第２次解线性方程组的迭代解法2。1例题Ｊacoｂｉ迭代法算法：流程图如下。高斯-塞德尔迭代算法：计算步骤与流程图与雅可比迭代法大致相同，只是一旦求出变元的某个新值后,就改用新值替代老值进行这一步剩下的计算。例2．1用Ｊaｃobi迭代法解以下方程：解：编制迭代计算的Ｍ文件程序如下:％Jacｏｂi迭代法求解例3．1%A为方程组的增广矩阵clｃ；A=[１0-2-13；-2１0—11５;-１-25１0]MAXTＩMＥ＝50；％最多进行50次迭代eps=１e-５;％迭代误差[n，m]=ｓize（A);x=zeroｓ(n,1）；%迭代初值y=ｚｅrｏs(n,1)；k=0;％进入迭代计算disp(’迭代过程X的值情况如下：＇)disp（’Ｘ='）;while１dｉｓp(x＇);ｆｏri=1:１:ns=０.0；forj=1：1：nifj～=ｉs=s+A（i，j）*x（j）；eｎdy(i）=（A(i,n+1）—s）/A（i,i)；endｅｎdfori＝1:１：nｍaxeps=max（0，ａbs(x(i）-ｙ（i)））；％检查是否满足迭代精度要求eｎｄifmaxeps＜＝ｅps％小于迭代精度退出迭代foｒi=１:1:nx(i）＝y（i）；％将结果赋给xendreturn；ｅndfｏri＝１:1:n％若不满足迭代精度要求继续进行迭代x（i）=y（i）;ｙ（i）=０．0；ｅｎｄk＝k＋1；ｉfk〉MAXTIＭＥ％超过最大迭代次数退出error(’超过最大迭代次数，退出＇)；ｒeturｎ；enｄend运行该程序结果如下：A=1０—2－13－２１0－115—１-２510迭代过程Ｘ的值情况如下：X=0000．3００01.50002。00000．８000１。76０02。660０0。9１801。92６02。86４0０。９7161。９7０02．９5400．98941。98972。98２30。99621.99６１2。９9３80。99８6１。9９8６2.9977０.99951。99９52。９9920。999８1.９9982．99９70．99９91.9９99２．99991.00０02．0000３。00０01.０0００2.00003.０000容易看出迭代计算最后结果为：.例２。２用Gａｕss-Seｉdel迭代法计算例２。１并作比较.解：编制求解程序Gauｓs_Seｉdeｌ。m如下：％Gaｕss_Seidel。m％Gaｕss_Ｓeidｅl迭代法求解例2。2％A为方程组的增广矩阵clc;ｆorｍatlonｇ；A＝［１0—2-13;－210-115；—1-25１0］[n,m]=sｉze（A);％最多进行5０次迭代Maxtime＝50；％控制误差Eｐs=1０E—5；％初始迭代值x＝zｅｒoｓ（1，n)；disp(’x=’）；％迭代次数小于最大迭代次数,进入迭代forｋ=1：Maxtiｍｅｄiｓp（x）;ｆori＝1：ns=０.0;foｒj=1:nifi～=js=ｓ＋Ａ(i，j）＊ｘ(j);％计算和ｅnｄenｄｘ（i）=（A(i，ｎ+1）—s）/A(i，ｉ)；%求出此时迭代的值ｅnd%因为方程的精确解为整数，所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件ｉfｓum（(x—floｏr(x）)。^２)＜Epsｂrｅak；ｅnd；enｄ；Ｘ=x；disｐ('迭代结果：')；Xｆoｒmａtｓｈort；完成后直接在Matlab命令窗口中输入Gaｕss＿Seidｅｌ回车后可得到如下结果:〉>Gausｓ_SeiｄelＡ＝10—２-1３—210-115-１—２51０x=0００0．３0001.５６00２。68400．8８001.９4448０0000０0０002.95387２0000000000。9８4283２0００00000１．99224３84００00００02。９93754１７600000０0。９9782４1８5600０００１．9989402547２0000２。99914０9390080000.９９971.99985４５22８69760２．99988２2381168６４０．9999５9128385638１。9９99800４94８8８１42。９9９９８384５４726５40。９99９9439４4450281.99999７2６3436２712。9999977842６３５140．99999923１１1360６１.99９99９６2464９0722。９9999969608２3５０0。９999９98945380４91.999９9994851584５2。999９9９9５8313948０。9９99999855345６4１．９99999992938３０８2.999999９942８223６0．99999９998015885１．９9999999９０3１401２．9９999999921５7３7０。9９999999９７27８5４1。9９９9999998６71452。９999９９999８９2４2９0。99999９９9９96２6721．99９9９9９９９9817７８2。9９999９99９9852460。99９9999999９4880１.999９99９999975０12。999999999９9797６0．99999９99９999２981.9999999999996５72.９9９99999９9997２２0。９99９99999999９04１。99999999999９95３2.999９99９999999620．９99999９999９99871。9９99999999９99９4２.9999９９９9９９9９9９5０.99999９999999998１.99９9999999９9999２．99９9999999999991。00002。0００03．０００0迭代结果：Ｘ=1２３可见对此题Ｇaｕss＿Seidel法的收敛速度还是很快的.２。2Matlab迭代解法常用函数介绍１．计算范数命令ｎormnoｒm(Ａ，option）当A为向量V时:nｏrm（Ｖ，1）为A的1—范数.noｒm（Ｖ，2)或noｒm（Ａ)为A的1—范数。nｏtm(V,Ｐ)＝sum(abs(V）.＾P）^(１／P）．norm（V，iｎf）=ｍaｘ（abs（V））为Ａ的无穷范数．norｍ(V，—inｆ)=min(abs(Ｖ））为A的最小范数．当Ａ为矩阵A时：noｒｍ（Ａ,１）为Ａ的列和范数。norｍ（A，２）或norm（A）为A的谱范数．norm(A，ｉnf)为A的行和范数。norm（A，'fｒo'）时为A的F—范数.2.计算矩阵的谱半径r＝ｍax(aｂs（eiｇ（A))）2.3实验作业实验名称：解线性方程组的迭代法．实验目的:1.掌握解线性方程组的雅可比迭代和高斯—塞德尔迭代算法；２．培养编程与上机调试能力．实验内容：应用雅可比迭代和高斯—塞德尔迭代算法解线性方程组：实验要求：(1)用Maｔｌaｂ语言编写雅可比（Ｊacoｂｉ）迭代法程序,并且选择不同的迭代次数，观察输出结果;(2）用Ｍａtlaｂ语言编写高斯-塞德尔(Gａuss—Sｅidel)迭代法程序，并且选择不同的迭代次数,观察输出结果．ﻬ第３次非线性方程求根二分法：算法流程图迭代法：算法流程图牛顿迭代法:算法流程图3。１例题例3。1判别方程的实根存在区间，要求区间长度不大于1，并求出最小正根的近似值,精度.解:先用画图的方法来粗略估计其根的范围.＞>ezploｔ(＇x^3-3*ｘ+1'）；%亦可用fplot命令＞〉griｄｏn；观察图可知其根分布区间大概分布在(－2，－1）,（０,1)和(１,2）中．编写二分法求解最小正根的近似值程序如下:foｒmａtloｎg;f=ｉnｌｉne（＇x＾3-3*x+１’）；a＝0;b＝1;Ｅps＝1Ｅ－5；forｋ=1:50A（k)=a;B(k）=b；ya=ｆeval(f,a）;yb＝fevａｌ(f，b）;tｅmp=(a+ｂ)/2；X（k)=teｍｐ；ｙt=ｆeｖaｌ（f，temp）；F(k）=ｙｔ;ifａbs(ｙt）〈Epｓｂｒｅak;endｉfyｔ＊ya〈0ａ=a；b=ｔemｐ；eｌseiｆｙt＊ｙb<0a=temｐ;b＝ｂ；eｎdeｎd;diｓp（'ｋａ（ｋ)ｂ(ｋ）x(k）ｆ(ｘ)')；H=［[1：k］',A’，Ｂ'，Ｘ＇，F']；disp（Ｈ）；ｄｉsｐ（'ｘ＝'）；diｓp（X（k))；diｓp（'y='）;disp（yt)；fｏrmatｓhort以下为运行结果：ka（ｋ）b(k)x（k)101.00000。50００2００.５００００。2５0０30.２５0０0.５０000。3750４０。2500０。3７５0０.３12５50。3１２50。375０0．343７50000０0000060．３4３75000０0００0000．375００．35937500００000007０。34３７500０000０00００。3５９３７500０0000000．３50０８0.３437５00000000000.35000.３４76５6２500000009０。3437500000００0000。347656250000０0０0.３457100．３４570．３4７６5625０0０00000.34６6７96８750000０110.３４６67９６８750０0０00．3４765６２5000０0000。347167968750000１20。34７１６79687500０００．347656２500０00０00。347４1210９３７50０01３0．34７16796８7500000.34741２1093７50000．347２90０３9０6２5001４0。347２900390６250０0。3４741210９３７５００0０。347３51０7４21875015０.34729003906250０0.３473510742１87５00.3４73205５6６406２51６0。34７29003906２5000．34732０556６406２５０.34731７０。347２9０0３9０6２5000.３４730。3４7２9766845703１ｆ（ｘ)—0.3７500．2656２500０００0000-0。0０0００。5000０．43７5-0。914１—0.874０－0.７3380。5２２2０.0４000。2650—0.44７00。9015—０.９67８－０．5734－０.０949－０．1613x＝0．34729７66845703１y=—3.464２2161３125３1８e－0０6从结果中分析可得到二分法求解该方程的根为x=０。３47297668４５7031.例3.2采用不同的迭代法,求方程在（１,２)内的近似根．解：构造迭代格式如下:(1）对应求解程序程序编制如下:f=inｌine（’x－x＾３-4＊x^2+1０'）％作内联函数ｄisp（＇x=');ｘ=feval(ｆ，1.5)；％迭代初始值diｓp(ｘ);ｉ=1；ｗｈilｅ1％进入迭代i=i+1；x=feｖaｌ(f,ｘ);ｄisp(x)；ifx>1E10;％此为x发散bｒeak；ｅnd;enｄ；ｉ，x运行结果如下:ｆ=Inｌｉnefｕｎｃtion:f(x）＝ｘ—x＾３-4*x^２＋1０x=－0。８７５06．７３24-469。７200１.０275e+０0８—1.０84９e+0241.2７71ｅ+07２i=6x=１.2771e+07２迭代6次后ｘ的值大得令人吃惊，表明构造的式子并不收敛．（２）：求解程序如下：f＝ｉnline（’(10／x－4＊x）^（1/2)’)dｉｓｐ（’x=’）；x=feval（f,1。5）;diｓｐ(ｘ）;Eps=1E-5;ｉ＝1；while1x0=x；i=i+1；x=ｆeｖal（f，x）;disp(x）；if~iｓrｅal(x）%不是实数不进行迭代breａk;ｅndifx＞1Ｅ10;%发散不进行迭代ｂreaｋ；end;ifａbｓ(ｘ-x0）〈Epsbｒeak；eｎdenｄ;i,ｘ运行结果如下：f＝Inlinefｕｎction：ｆ（x)＝（10／x-4＊x)^（1／2）ｘ=０。8１652。9９690。0００0+2。9412ｉｉ＝3ｘ=0.0000+２。9412ｉ由以上结果可知当进入第三次迭代时，出现了复数,这并不是我们期望的，因此该式并不收敛．(３)程序如下，这里将精度提高至1E－10：fｏrmatlong;f=inlｉne（＇1／2*(10-x^3）＾（1/2)＇)diｓp（'x＝'）；x＝fevaｌ（f,1.5)；disp(x);Eps=1E—10;i＝１；while１x０=x；i=i+1;x＝fevａl(ｆ，x）;ｄiｓp（x)；ｉｆ～iｓreal(x）ｂreａk；enｄｉfx>1E1０；ｂreaｋ；eｎd；ifabs(x-x0）〈Epｓbreak；endenｄ;i,xｆorｍaｔshｏrt；运行上述程序得到如下结果:f=Ｉｎlineｆunｃtion:f（x）＝1／2＊（10－x^3)^(１/2）x=１.2８6９537676233７51。４578１.3454583７40２32941.３75１7０2５２8160381.3６33１。36784696７5921331。363８870038８402１1.３6５9167３3３900４01.３6４878217193６771.3６５41０0611６９9５７1．36５31。3６52772０852４4791。36５21．3652423８37188３9１．365223680225２8２1。3６5２３３25５742５０0１。365228353４626２71．36523086３２4３６3７1。36522９57８3339591.3652302361５８１81１．3６5２2９89９３7７7331.36５23０07179６2911.3652299８3524６741.3６52３00２8716３23１。３652300０55７995０1。３6５230０１７4248771。３65230０1１36０7331.3６５2３０0１4４653４0１．36５２3００１28759０1１.3６52300136８963２1.36５2３00132７303４1.36523001348６３161.3６５２30０1３37７1241。365２300１３4３3026ｉ=３4x=１.36523００1３43３026从结果可知该式子收敛但速度较慢.（４)程序将上例中的内联函数改为f=inline（'（10/（4＋x)）^(1／２）’)即可,运行结果如下:f＝Inｌinefuｎctioｎ：f(ｘ）=（10／（4+x)）＾（1/２）ｘ＝1.3４8３９９724９2６４841.３67３7６3７19９1２83１．364９5701５４024８7１.3６５２64748113４421.3652２559416０52５1．3652305７567３4341．365２2994１８78183１。36523002２515５691.3652３001２2561２2１。3６5230０135614251。3６52３001３395３521．３65２300134１6482i=12x=1。3６523０013４1６482从结果中可知该式子收敛很快，只需12次迭代即可达到(3）的精度．（５）仍将（３)中的内联函数替换成f=iｎlinｅ(’x—（x^３+４＊x^2－１0）／(3*ｘ^2＋8*x）＇),运行后得到结果如下：f=Inlineｆuｎctioｎ：f（x）=x-（x^3+4*x^２—１0)／（3*x＾2+8*ｘ)x=1。373３33333３33３331。3652620１48７46２71。３６52３００１3９１61471.36523001３41409７1。3６52３00１３414097i＝5x＝1．36５23０0134１40９7该式子收敛相对(4)式更快，结果也很准确.例3．4用Neｗton法求在0。7附近的根,计算结果保留６位有效数字。解：构造符号函数：〉〉x=ｓym（'x')〉>f=sym(＇ｘ＊sin（x)-0.５＇）求符号函数导数>>ｄf＝diff（f，ｘ）构造Neｗｔon迭代公式>＞ＦX=x-f/ｄf%以下为输入结果x=xｆ=x＊sin(x）－0．5dｆ=siｎ（x）+x＊cos（ｘ）ＦX＝x－(x＊sｉｎ（x）-．5）/(sin（x)+x*cos(ｘ))进入迭代计算，计算前10次迭代值：fｏrmａtlong;Fx=inｌinｅ（FX)；x0=０.７；forｉ=１：１0disｐ（x0);x0=feval（Fx,x0)；eｎdｘ0fｏrmａtshｏｒt;输出ｘ的中间结果为：0。70０00.7４830。7４560。７4１０0。74910.７4910.74910．７４910.７491０.74９１x0=0.7491可见法迭代5次即趋于稳定．3。2Matlaｂ非线性方程求根的命令1。代数方程组的求根ｒootｓr＝rｏｏts(P）P为代数方程的系数向量,从高次到低次排列,该命令只能求出一个一元多项式的根，返回所有复数和实数根.2。求零点fzerox=fｚeｒｏ（F,x０,oｐtｉｏｎ）F为函数和字符表达式、内联函数或M文件的函数形式，x0为零点的大概位置,optｉｏｎ为可选项,可参阅其他资料得到详细说明。使用该命令前,前配合plｏt、ezpｌot和fpｌot先画出函数图形,再猜测ｘ0的大概值。3求方程组数值解的命令fsolｖex＝fｓolｖｅ(’fun’,x0，optionｓ）fun为M函数文件名，,它是探索方程的起点,输出为与ｘ0同维的向量或矩阵,是方程组的数值解。对于更详细的设置opｔion,可用ｈeｌpfsolve获得详细说明。4.非线性方程组的解析解ｓolvesolｖe('eqn1',’eqn2',…,'eqnN＇)对Ｎ个方程采用默认变量求解.soｌve（＇eqn1’，ｅqn2’,…，＇eqnN',’vａr1,ｖａr２,.。，varN＇)对N个方程的var1，var2，．．，varN变量求解．Ｓ=ｓｏlｖe('eqｎ1＇，eｑn2',…，'ｅｑnＮ＇，’ｖar1＇，'var2'，。．,＇ｖarＮ’)对Ｎ个方程的vａr１,vａr2,．.,vaｒN变量求解.［x1,x2，…，ｘn］=ｓｏlve('eqn1’，eqn2’，…，＇ｅｑnN＇,'vａr1’,＇var2＇,.．，’vaｒN'）将求解分别赋给变量x１,ｘ２,…,xｎ。3.３实验作业实验名称:解非线性方程的迭代法.实验目的：掌握解非线性方程的迭代法；实验内容:求方程的根。实验要求:用牛顿迭代，同样计算到为止．输出迭代初值及各次迭代值和迭代次数．第４次插值法拉格朗日插值:算法流程图4。1例题例4．1已知，用线性插值法求的近似值．解：Matlab中有直接进行线性插值计算的命令interp１，直接使用intｅrp１命令即可.〉＞x=［４9];>〉y=［23];>〉ｆ=inｔeｒp1(x,y,7，’ｌinｅar'）％选项使用线性插值ｆ=2.６00０故插值计算结果为例４。2,用二次Laｇrange插值求的近似值。解:求解过程描述如下:formatloｎg；％显示15位%输入初始数据x0＝［4916];ｙ0＝［2３4];ｘ=7;％插值点n=lｅｎgth（ｘ0)；%取长度％初始计算s=０;％进入公式计算forj=0：(ｎ-１）ｔ=1;fori=0:(n—1）ifi~=jt=t*（x—x0（ｉ＋１））/（x0(j+1)—x0(ｉ＋１）);ｅndends=s＋t*ｙ０(j+1）;eｎｄs％显示输出结果forｍatshoｒt;程序输出结果：s＝例4.3设有函数在［－5，5］上取等距节点上进行Lagranｇe插值.解:采用n次进行插值，为使计算结果可视化，画图显示插值结果,编制如下程序:％输入初始数据clc；x0＝—5:5；ｙ0=1./（１+x0.^2)；x=－5：0。0５：5；％插值计算点m=length(x）;ｆork＝1:ｍ%分别对每一点进行插值tｘ＝x（k);%插值点n＝ｌength(x0）;s＝０;%进入迭代计算过程forj=0:(n－1)t＝１；fori＝0：（n－1）ｉfi～=jt=t＊（ｔｘ-x0（ｉ+1））/(x0(ｊ+1）－ｘ0(i+1））；ｅndｅnｄs=s＋t＊y0(j+1);endyｘ（k）=s;end％画图显示结果ｆigure；ｐlot(ｘ,yx,＇：m＇,＇ＬineWiｄth’,2）;holdon;pｌｏt（x０,y0,’-r','ＬiｎｅＷｉｄtｈ',0。５）;gridoｎ;tiｔle（＇Laｇranｇe插值时的Ｒuｎge振荡现象'）;xlaｂｅl（＇x轴’）;ylabel（’y轴’）；lｅgend(’插值曲线’,'原曲线')；运行程序输出结果如下：由图可知确实在附近发生了振荡,误差很大。例4。4.设有函数在[—5,5］上取等距节点上进行分段线性插值并且绘制图像。Matlab程序如下Funcｔiony=dｉｖ_ｌineaｒ（x0，y0，x，n）ｆorｊ=1:lenｇｔh（x）；fｏri=1:n-1ｉf(x〉＝x0（ｉ）)&＆(x〈=ｘ０（i+１）)y=(x—ｘ0（ｉ+1）)／（ｘ0(i）-ｘ0(ｉ+1））*y0（i)+（x－ｘ0（i））/(x０（i+1）—x０(ｉ)）＊y0(ｉ+1）；elsecontｉnuｅ；end%endcoｎdｉtｉｏnalstａｔemｅｎtｅndｆuｎｃtiｏny=ｆ（x）%原来的函数ｙ=５。/（x。^2＋1）；x0＝—5：1:５；%插值节点w=length（x0)；％插值区间长度n=w－1；％节点个数forｘ=－５：０.01:5%加密了的插值节点y=dｉv_lｉneaｒ（x0，f（x0），ｘ，ｎ）；％分段插值后的点的纵坐标hｏldｏn;plot（x，ｙ，‘r’);plot(x,f（ｘ），‘b');ｅnd运行节点画图如下例4.5设有函数在[－5,5］上取等距节点上进行三次样条插值,并且绘制图像。注解：我没有理解此程序。上实验课程的时候可以选择使用Matlaｂ已经封装了的函数y=spｌiｎｅ(x0,y０，ｘ）直接进行计算并且画图.解：％X为１阶导数横坐标向量％Ｙ为1阶导数纵坐标向量％dx0和dｘn导数边界条件fuｎctiｏnS=cｓfiｔ(X，Y，dｘ0,dxn）%3次样条插值函数？N=lｅnｇth(X)—1;Ｈ=ｄiff（Ｘ）;Ｄ=diff(Y)./H；A=H(2：N-1）；％边界约束Ｂ=２＊(H(1:N-１)＋H(2:N））C＝H（2：N);U=6＊diff（D）；B（1)=B（１)-H(1）/2；U（1)=U(1)-3＊(D(1）);B(N—1）=B（N-1）—Ｈ(N）/２；U(N—1）＝U(N—1）—３*(－D(N））;fork=2:N-1temp=Ａ（k-１）/B(k—1)；B(k)=B（k）－ｔemp*C（k－1）；U(k）=Ｕ(k)－teｍp＊Ｕ(ｋ—１);endM(N）=U（N－１）/B（Ｎ－1）;fork＝Ｎ-2:—1：1M（k+1）=(Ｕ（k）－C（k）＊M（k＋２)）／B（k）；endＭ(1）=3＊(Ｄ（1)-dx0)/H(1）—Ｍ（２）/2；M（Ｎ+1）＝3＊(ｄｘn—D（Ｎ))／H(N）-Ｍ（Ｎ）／2；foｒｋ＝0：N—１Ｓ(k+1,1)=(M(k+2)-M（k+1)）/6*Ｈ(k+1)；S（ｋ+1，２）=M(k+1)/2;S（ｋ+1,３)=D（k+1)-H（k＋1)＊(２＊M（ｋ+１）+M（k+2））/6；S（ｋ＋1,4)＝Y(k+１);ｅndcｌearallclｃｃｌｆx=－5：1:5；y＝５.／(x。＾2+1）;X=－5：1:5;Ｙ=ｙ;dｘ０=０。１4２；ｄxn=-0.142；S=cｓｆｉｔ(X,Y,ｄx0,ｄxn);ｘ1＝－5:0。0１:-4;y1＝ｐｏｌｙvaｌ（S（1,：)，x1－X(1）);x２=-4：0。01：—3；y2＝ｐolyval(S（２,：),x２－Ｘ(2))；x3=－3：0.01：-２;ｙ3=pｏlyｖal（S（３,：)，x3－Ｘ(3)）;x４=—2:0．01:-1；y4=pｏlyval(S(4,：），ｘ４－X(4）);x5=-1：０。0１:0;y5=polyvａl(Ｓ(5，：),x5—X（5）)；x6=０：０.０1：１;y６=polyval（Ｓ(6,：)，x6—X（6）)；x７=１：0.01:２;ｙ7＝pｏlyval（Ｓ(7，:）,ｘ7－X（７）);ｘ８＝2:0。０1:3;ｙ8=pｏlyval（Ｓ(8,：），x8-X(８）);x9=３：０。01:4;y９=ｐolｙvaｌ（S(９,：),x9-X（９)）；x1０=4：0。01:5;ｙ10＝polyval(Ｓ（10,：)，x10－X(10)）；ｓubplot（２，１,１)；pｌot（x１,ｙ1，x２，ｙ2，x3,y3，ｘ4,y4，x5,y5，x６，y6,ｘ7，y7，x８，ｙ8，x９,ｙ9，x１０,y10，X,Y,‘＊'）ｇrｉdonlegend（‘样条插值');x=-5：１：５；y=1．/(x．^2+1);subplot(２，1,2）;ploｔ(ｘ，y,‘r’）gridｏnlegend（‘原来的函数’);４。2Matlab插值函数介绍1.一元函数插值ｉｎteｒｐ1y=iｎterp1（x0,y0,x，'ｍethoｄ’)ｘ0和y0为已知的两个同维向量，ｘ为插值点,可为向量和数值或者矩阵,'ｍethoｄ'含义如下:’nｅａreｓt’——最近插值,即用直角折线连接各样本点．’linear'——线性插值,依次边接各样本点，为默认选项。’pchip＇（或’ｃubic’）——分段一次多项式Heｒmite插值曲线,具有一阶导数连续性．’splｉne＇-—三次样条插值．２。三次插值ｐｃhiｐy=pchip（x0，y0，ｘ)与ｙ=ｉnterp1(x０，y0，ｘ,’pｃｈｉp’)等价。3．三次样条插值ｓpliｎｅy=splinｅ（x0，y０,x）它与y=ｉnｔerｐ1（ｘ0，y0，ｘ，＇sｐｌinｅ’)等价．4.二维插值inteｒｐ2ｚ=interp2（ｘ0，y０，ｚ0，x,y，＇metｈｏd'）（x0,ｙ0，z０）为原始数据点，x和ｙ为独立变量，＇meｔhod’的参数与ｉntｅｒp１中的相同．本函数主要用于曲面插值．５．其它插值函数三维插值interp3:三维插值函数．interpft()：利用FFT（快速傅里叶变换）的一维插值函数。inｔｅrpｎ(）：高维插值函数．cｓape（)：可输入边界条件进行插值．４.3实验作业实验名称:插值实验目的：1．学会利用Mａtlａb实现拉格朗日插值；２。学会利用Ｍaｔlab实现分段线性插值；3．学会比较拉格朗日插值与分段线性插值．实验内容1：已知函数在下列各点的值为试用次拉格朗日插值插值多项式及分段线性插值函数对数据进行插值．用图给出，及．实验内容2：设有函数将［-5,5］４等分，６等分，8等分，1０等分（即取等距节点)进行三次样条插值（直接使用4.2中所列的Matlab函数），并且绘制图像。解释插值结果

第5次利用最小2乘法进行曲线拟合5．1例题例5。1设有如下数据:1345678９１01０5４2１1234利用最小二乘法求该组数据的多项式拟合曲线．解：先输入数据并画图判断曲线的大致形状：〉〉x0=1:１0;ｘ0（2）=［]x0=1３45６78910＞〉y0=[1０５42１1234］y0=10542１1２34>>fｉguｒe,pｌot（ｘ０，y0,＇—r’,’LinｅWiｄth＇，3）,gridon，tｉtle（’数据观察图＇）；holdon,ｐlot(x０,y0，'ｏｇ＇）；axis（［011011］)；观察到图形大概为一条抛物线,因此它的最高次数为2.构造其正规方程的程序如下:ｘ0=１：10；x0（２）=[］;%初始数据％x为数据点的横坐标，y为数据点的纵坐标x=x0;y0=［1０５421123４］；y＝y0;m=2;%最高次数为2ｎ=lenｇth(ｘ）；b=zｅroｓ(1,m+1）;f＝ｚeｒｏｓ(n，m+1）;％f为正规方程的系数,初始为0fork＝１：m＋1f(:，k)=x’.＾(k-1)；ｅnda=f'＊f;b=ｆ’＊y'；%解方程，得到多项式由高到低的系数所构成的向量cc=ａ\b;c＝flipud（c)；dｉsp（ｃ)；运行上述程序得到数据的系数如下（从高次到低次排列)：0.26７５70６646２9484-3．68１9１3。45９6638655460９8至此,得到了二次多项式。画出该2次多项式的图:〉>c＝c’；〉〉x=0：0．1：11；〉>Ｆ=c（1)*x．^2+c（2）＊x+c(3）;＞〉plot(ｘ0,y0，＇ｏｋ'，＇LiｎeWidth’，2）,griｄｏn;〉〉ｈoldoｎ；＞>ｐloｔ（x，Ｆ，＇—ｒ',’LineWidth＇，１）；＞＞tiｔle（’拟合后的效果图’）;〉〉aｘｉs(［0１101１］);拟合后的效果如下(带点的是原始数据,红色的是拟合２次曲线）：例5.2设如下数据:0０．25０。500．７51。001．00001.２８40１.6487２.11702.7183求其拟合多项式．解：按解题步骤,先求其一次多项式，由于以下很多拟合都要用到拟合系数，我们把例5。1中的主要部分抽取出来作为函数使用,以节省代码编写量.fｕncｔionｃ=mｙpｏlyfit(x,ｙ,m）％myｐolyｆit。m％ｘ为数据点的横坐标％y为数据点的纵坐标%最高次数为ｍ次%返回结果为c,c按次数从高至低排列成向量形式n=lｅｎgth（ｘ）；ｂ=zｅros（1,m＋1)；f=zeros（n，ｍ+1）；％f为正规方程的系数，初始为0ｆoｒk=1:ｍ+1f(:，k）=x'.^（ｋ-1）;eｎda=ｆ＇＊ｆ；b=f＇＊y'；%解方程，得到多项式由高到低的系数所构成的向量cｃ=a\b；c=flipud（ｃ）；％翻转另存为myｐolyfit.m后，以下进入拟合阶段:>〉x=0:０。25:1；%输入初始数据〉>y=[11．2８41。64872．1１702。７183］；>>m=１;%先作最高拟合次数为1>〉C=myｐolyｆit（x，y，m）％得到拟合多项式系数C=1.７0０10.８9９６８０00００00００0＞>ploｔ(x，y,'ｏk’,’LiｎｅＷiｄｔh’，2）;%画出拟合的效果图〉〉gridon;＞〉tｉtle(＇最小二乘一次拟合效果图＇)；〉>x０=0：0.1:1；〉>ｙ0=C(１）．＊x0+C（2）；>〉holdoｎ，ｐｌｏt(x0,y０，’－ｒ'）;最小二乘的一次拟合效果如下：＞>m=２；%改为拟合二次多项式>＞Ｃ=mypolyｆｉt(ｘ,y，m)%求出拟合系数C=0．8４3６57142857134０．８６4７１。０05２＞〉fｉgurｅ；％画图得到拟合的效果图＞〉ｐlｏt(x,ｙ,’oｋ＇，'LｉneＷiｄtｈ＇，２）;〉>grｉdｏｎ；＞＞ｔitle(’最小二乘二次拟合效果图')；＞>x０=０：０.1：1;>〉y0=Ｃ(１).＊x０。^２+Ｃ（2）.＊x０+Ｃ(3）；〉＞holdｏn，pｌot（ｘ0,y０，’－r＇);最小二乘的二次拟合效果如下：５．2Matlab数据拟合命令介绍1．数据的多项式曲线拟合pｏlyｆitp=polyｆiｔ(ｘ，y,m）x，y为样本点数据，m为拟合的最高次数,ｐ为拟合的多项式．求出拟合系数后，可配合poｌyｖaｌ（p,x０）计算多项式在x0的值，x0可为向量或数值.另外还可用polｙ2ｓtr（p）将表达式的形式显示出来。以下再介绍几个与多项式有关的命令．2.多项式创建ｐｏlyp=poly（A）将向量A转化为多项式．３．多项式转换poly2ｓｔrpｏｌy２str(p,'t＇)将多项式的数学形式恢复出来.４.多项式运算conv(p1，p2)返回p１和p2的乘积.［q，r］＝deｃｏnｖ(ｐ1，p2)q为ｐ1除以p２的商多项式系数，r为余数．ｋ=pｏｌydｅr（p）多项式求导．pｏｌyvａl(ｐ，x0）求多项式p在x0点的值，x０可为一个数值或向量。5。Ｍａtｌａｂ的非线性拟合命令isｑｃurvｅfｉｔ():该函数位于优化工具箱内．isqnoliｎ（):位于优化工具箱内．nｌinfit（）：位于统计工具箱内.5。3实验作业实验名称：最小二乘法实验目的：学会利用Mａtlab实现最小二乘法.实验内容:设有观察数据－2．0-１．5-1．0—0。５０0.511。52．09。125．５82．961。４7１.01．513．05．４９６9.19用最小二乘法求出合适的拟合曲线,并且画出图。第6次数值积分复化梯形公式：算法流程图复合辛普生公式：算法流程图变步长梯形公式：算法流程图6.１例题例6.1给定积分，分别用梯形公式、公式、公式作近似计算．解:先输入主要初始参数〉>a=０．5;〉＞b=1；>>f=inlinｅ('x^(１/2)’）；％梯形公式〉＞I1=（b－a）/２*（feval（f,a）+fｅｖal(ｆ,ｂ）)I1＝0．426776６95２96637%simpsｏn公式>〉I2＝(b—a）/6*（ｆｅval（f，a）+４*ｆeval（ｆ,(ａ+b）／２)＋feval(ｆ，b）)I2=０。432５%Coｔｅs公式(n=４)〉＞ｔc=０;＞〉C0＝[７321２327];>>ｆｏrｉ＝０：4ｔｃ＝ｔc+C0（i+１）*feｖａl（f,a+i＊(b—a）／４)；enｄ>>I3=（ｂ—a）/9０*tｃI３=0.4３７6％准确值>＞I=int(char(ｆ)，ａ，b)〉〉ｖpa(I）Ｉ=—1／6*2^(１/２）+2/3aｎs＝0.43５４596505例6．２对积分,为使其精度达到.若用复化梯形公式，应将［０，1］多少等分？若用复化公式，应将［０，１]多少等分？解：直接按余项计算即可。复化梯形公式的余项为::复化公式余项为:对于,在数学分析中，我们已证得以下不等式成立:．直接利用上述不等式关系解答本题.先输入误差精度：>〉Eps=1E－4Eps＝1．0000ｅ－０0４复化梯形公式>〉h1=ｓqrｔ（Eｐs／aｂs（—（1－0）/12*1／(2+1))）%先求出步长h1=0.０000>〉N1=ceiｌ（1/h1)％向上取整,得到等分区间数Ｎ１=1７故可将区间17等分即可达到所要求的精度.(2)复化公式〉>h2＝power(Eps/aｂs（-(1－０）/180＊1/（1＋4))，１/4）%先求出步长h2=０。5477225575051６６＞＞Ｎ２=ceｉl(1/h2）％向上取整,得到等分区间数Ｎ2＝2故可将区间2等分即可达到所要求的精度.扩展：Maｔｌab中复化梯形公式命令为I=trａpｚ（ｘ，y),复化公式命令为quａd()．Matlaｂ中有四个取整函数,分别为ceiｌ（),ｆloｏr（），fix(），round()，分别表示向正无穷大方向取整、向负无穷大方向取整、向靠近零方向舍入和四舍五入.例６。3对积分，利用变步长方法求其近似值，使其精度达到．解：利用变步长法前,先建立三种变步长复化积分公式的函数.注意在Ｍatｌab中直接用siｎ（0)／０得不到1,，因此解此题时,我们改用求极限的方法得到函数值，此函数名为limit(）。先建立三种复化公式的函数文件，它们分别为复化梯形公式trａp。m、复化公式为simpson。m、公式为cｏｔes.m,三个函数的源程序如下：复化梯形公式ｔｒap。mｆunctionT＝trａp(f,ａ,b,n)%复化梯形公式求积分值%ｆ为积分函数％［a,b］为积分区间%n是等分区间份数h=（b-a）/n;％步长T=0;foｒk=１：（n－1)x０=a+h*k;T=Ｔ＋limit（f，x0);endT=h*（ｌimit(f,a)+limit(f,b))／２+h＊T；T=doｕｂle(T）；（2）复化公式simｐsｏn。m：fｕｎctionS=simpsｏn（f，ａ,b，n)％Ｓｉｍｐson公式求积分值%f为积分函数％［a,b］为积分区间%n是等分区间份数h=（ｂ—a)/（2＊n);％步长s1=０；s２=0;forｋ=１:nx０=a+h*(２＊k—1）；s1=s1＋liｍit(f，x0）；eｎdforｋ=１:（n－1）ｘ０＝a＋h*2*k;ｓ2=s2＋ｌｉｍit(f，ｘ0）;enｄS=h＊(limｉｔ(f,a）+limit（ｆ，b)+4＊s1+2*s2）/3；S=doｕblｅ（S);(3)复化公式cotes。ｍ：functｉonC=cｏtｅ（f，ａ,ｂ，n)%复化cｏtｅｓ公式求积分值%f为积分函数％［ａ,b]为积分区间％n是等分区间份数h=（b—a）/n;%步长C＝0；ｆori=1:（n-1)x0=ａ＋i＊h；C＝C+14＊limit(f,x0）；endｆｏrk=0:(ｎ—１）ｘ０=a+ｈ＊k;s=32＊ｌi