数值分析上机第三次
昆明理工大学数值分析上机报告3
数值分析实验报告（数值积分）姓名：学号：2006231011专业：材料学学院：云南省新材料制备与加工重点实验室授课教师：昆明理工大学06工科硕士 《数值分析》上机实验报告专业： 材料物化 姓名： 学号： 2006231011 任课教师： 作业完成实验室：实验内容：1．题目/要求：1、 利用Lagrange 插值公式()k n ki i ik in k n y x x x x x L ⎪⎪⎪⎭⎫ ⎝⎛--∑=∏≠==00 编写出插值多项式程序； 2、 给出插值多项式或分段三次插值多项式的表达式；3、 根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何； 对此插值问题用Newton 插值多项式其结果如何2．作业环境(包括选用的程序语言、运行环境) Visual C++ 6.03．数学（理论背景）描述1. Lagrange 插值多项式定义：若n 次多项式上满足条件个节点在,1),1,0)((10n j x x x n n j x l ⋅⋅⋅<+⋅⋅⋅=： 当0)(1)(=≠==k j k j x l j k x l j k 时；当时， 其中n k j ⋅⋅⋅=1,0, ，就称这1+n个次多项式)(),(),(10x l x l x l n ⋅⋅⋅为节点n x x x ⋅⋅⋅10,上的n 次插值基函数。
插值多项式可 表示为：)()(0x ly x L knk k n ∑==，称为Lagrange 插值多项式。
2.分段线形插值就是通过插值点用折线段连接起来逼近).(x f 设已知节点b x x x a n =<⋅⋅⋅<=10上的函数值n f f f ⋅⋅⋅,,21，记k kk k k h h x x h ma x ,1=-=+求一折线函数)(x I h 满足：[][]为分段线形插值函数上是线形函数，则称在每个区间记)(,)(3),,2,1,0()(2,,)(11000x I x x x I k f x I b a C x I h k k h k h h +==∈4．数值计算公式Lagrange 插值多项式()k nki i ik in k n y x x x x x L ⎪⎪⎪⎭⎫ ⎝⎛--∑=∏≠==00； 分段线形插函数，在每个小区间[]1,+k k x x 上可表示为：1111)(++++--+--=k kk kk k k k h f x x x x f x x x x h I （1+≤≤k k x x x ）在整个区间[]b a ,上为：)()(0x lf x I jnj j h ∑==其中基函数),1,0,()(n k j x l jk k j ⋅⋅⋅==δ，其形式是：[][],.,,,,0);(,);0(,)(11111111⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧∉∈=≤≤--=≤≤--=+-+++---j j j j j j j j j j j j j x x x b a x n j x x x x x x x j x x x x x x x x l5．算法程序流程图Lagrange插值算法程序流程图分段低次插值算法程序流程图6．程序结构(程序中的函数调用关系图6．实验数据和实验结果（打印或用屏幕图形拷屏表示，可加为附页）（1）Lagrange 插值算法,62573238.0)596.0(=f 05422977.0)99.0(=f 分段低次插值算法,63270600.0)596.0(=f 10152800.1)99.0(=f （2Lagrange 插值算法16476189.0)8.1(=f 00126583.0)15.6(=f分段低次插值算法18160000.0)8.1(=f 00185000.0)15.6(=f7．讨论（包括题目要求的讨论和方法的适用性讨论）对于插值多项式L n (x)，当∞→n 时，L n (x)不一定收敛到f(x)，此时需用分段线性插值比L n (x)逼近f(x)好得多，题（2）就是一例。
数值分析第三次作业
26.解：(1).J 法：J ∴法收敛.GS 法：()11102210221101110122100210G B D L U ----⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=-=-=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦⎣⎦022023002-⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦()()()212322det 0232020,2,21G G I B B λλλλλλλλλρ--=-=--∴===∴=GS ∴法不收敛.()2.J 法：()()131231*********()12202101121101212012125det 412125550,,,1222J J J B D L U I B i i B λλλλλλλλλρ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+=--=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦--==+--∴===-∴=J ∴法不收敛.GS 法：()()()1312310220221101101122022022det 11002201J J J B D L U I B B λλλλλλλλρ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+=--=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥----⎣⎦⎣⎦⎣⎦--===⇒===∴=()1120111200011220212120021120014120G B D L U ----⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=-=-=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎣⎦01212012120012-⎡⎤⎢⎥=--⎢⎥⎢⎥-⎣⎦()()()21231212det 01212120120,12,121G G I B B λλλλλλλλλρ--=+=+=+∴===-=GS ∴法收敛27.解：()1010911102,702106A b -⎛⎫⎛⎫ ⎪ ⎪=--= ⎪ ⎪⎪ ⎪-⎝⎭⎝⎭A 为严格对角占优阵，J ∴法和GS 法均收敛.J 法的分量形式为：()()()11111,1,2,,i nk k k ii ij j ij j j j i ii x b a x a x i n a -+==+⎛⎫=--= ⎪⎝⎭∑∑J ∴法的迭代格式为：(1)()12(1)()()213(1)()321(9)101(72)101(62)10k k k k k k k x x x x x x x +++⎧=+⎪⎪⎪=++⎨⎪⎪=+⎪⎩取初值(0)0x =，J 法的数值结果是：迭代次数k()1k x ()2k x ()3kx 1 0.900000 0.700000 0.600000 2 0.970000 0.910000 0.740000 3 0.991000 0.945000 0.782000 4 0.994000 0.955500 0.789000 50.9955500.9572500.791100GS 法的分量形式为：()()()111111,1,2,,i nk k k ii ij j ij j j j i ii x b a x a x i n a -++==+⎛⎫=--= ⎪⎝⎭∑∑GS ∴法的迭代格式为：(1)()12(1)(1)()213(1)(1)321(9)101(72)101(62)10k k k k k k k x x x x x x x +++++⎧=+⎪⎪⎪=++⎨⎪⎪=+⎪⎩取初值(0)0x =，GS 法的数值结果是： 迭代次数k ()1k x()2k x()3kx10.900000 0.790000 0.758000 2 0.979000 0.949500 0.789900 3 0.994950 0.957475 0.791495 4 0.9957475 0.9578738 0.7915748 50.9957874 0.95789370.7915787()123210,99,950A∆=∆=∆=∴对称正定，()1110000101001011100102010105020110000105J B D L U -⎛⎫⎛⎫ ⎪ ⎪⎛⎫ ⎪⎪ ⎪⎪ ⎪=+=-= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭ ⎪ ⎪ ⎪⎪⎝⎭⎝⎭212,310101111det()()00,10520201051()20J J I B B λλλλλλλλρ-∴-=--=-=⇒==±-∴=∴SOR 法的最优松弛因子为：[]2221.01282111/2011()opt J B ωρ==≈+-+-()10.01282opt opt L ρωω=-=对应的渐近收敛率为：R=-ln ()() 4.35654opt opt R L L ωρω=SOR 法的分量形式为：()()()()()111111,1,2,,i nk k k k iii ij j ijjj j i ii x x b a x a xi n a ωω-++==+⎛⎫=-+--= ⎪⎝⎭∑∑∴SOR 法（ω取最佳松弛因子）的迭代格式为：(1)()()112(1)()(1)()2213(1)()(1)3321.012820.01282(9)101.012820.01282(72)101.012820.01282(62)10k k k k k k k k k k x x x x x x x x x x +++++⎧=-++⎪⎪⎪=-+++⎨⎪⎪=-++⎪⎩取初值(0)0x =，SOR 法的数值结果是： 迭代次数k ()1k x()2k x()3kx10.911538 0.801296 0.770006 2 0.981009 0.954035 0.791074 3 0.995588 0.957822 0.791571 4 0.995785 0.957894 0.791579 50.9957890.9578950.79157928.一定收敛.证明：对于11122122a a A a a ⎛⎫=⎪⎝⎭，A 对称正定，()212211122122111221201,2,,det()0,iia i a a A a a a a a a a ∴===-对于J 法：121111121121222221000()01000J a aa a B D L U a a a a -⎛⎫⎛⎫-⎪ ⎪-⎛⎫⎪ ⎪=+== ⎪ ⎪ ⎪-⎝⎭- ⎪ ⎪⎝⎭⎝⎭122211212121,2121122112222det()0J a a a a I B a a a a a a λλλλλ-==-=⇒=22121122121122()1J a a a a B a a ρ∴=∴J 法收敛. GS 法：12221112121221122112212112222121122121122det()00,0()1G G a a a a I B a a a a a a a a a a a B a a λλλλλλλρ⎛⎫-==-=⇒==⎪⎝⎭-∴=∴GS 法收敛.∴对于系数矩阵对称正定的2阶线性方程组，J 法和GS 法一定收敛. 30.证明：由线性代数知识知：∃可逆矩阵使121s J J p B P J J -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦其中,i in n iJ R ⨯∈对应于特征值()121,2,,i s i s n n n nλ=++=()0B ρ=∴B 的所有特征值为0，120101,1,2,,10i i r r i s J R i s r r r n⨯⎛⎫⎪⎪⎪∴=∈=++= ⎪ ⎪ ⎪⎝⎭1i r =时，11i J R ⨯∈1i r 时， 0,i i r r i i J J R ⨯=∈12kkkk s J J J J ⎡⎤⎢⎥⎢⎥∴=⎢⎥⎢⎥⎢⎥⎣⎦最多迭代到第n 次，即k=n 时，10,0k k k J B PJ P -=== 设x *是Ax b =的精确解，误差向量()()k k e x x *=-()()()()()()110k k k k ke x x B x x B e B e--**=-=-===所以最多迭代到第n 次时，()()()00,k k k e B e x x *===所以结论成立31．解：（1）根据迭代公式(1)()()()k k k x x Ax b α+=--，有： (1)()()k k xI A x b αα+=-+ ∴迭代矩阵13212B I A ααααα--⎛⎫=-= ⎪--⎝⎭ 12132det()(1)(14)0121,14I B λααλλαλααλαλαλα-+∴-==-+-+=-+∴=-=-当{}()max 1,141B ραα=--时，迭代收敛111110121411141ααααα⎧---⎧⎪⇒⇒⎨⎨---⎪⎩⎩012α∴时，此迭代方法收敛{}()()m a x 1,141,00.441,0.40.5B B ρααααραα=---⎧∴=⎨-≤⎩ 0.4α∴=时，()Bρ最小，迭代收敛最快()12,n λλ为A 的特征值，11,1n αλαλ∴--为I A α-的特征值{}1()m a x 1,1n I A ρααλαλ∴-=--必要性：迭代收敛()1I A ρα⇒-111110211nαλαλαλ-⎧-⎪∴⇒⎨-⎪⎩所以必要性成立 充分性：()1111102022,1,2,11,1,2,()max 11i i ini i i ni nI A αλαλλλαλρααλ--=∴≤=∴-=∴-=-所以此迭代法收敛，充分性成立 （3） 1102αλ-时，111121,0()21,2n n in I A αλαλλρααλαλλλ-⎧-≤⎪+⎪-=⎨⎪-⎪+⎩根据图像，12nαλλ=+时，()I A ρα-最小33.解：()()()()()()()()()()()()()()()1()1121()111211111(1)()11111,k k k k k k x D L Ux D L b x L D U x L b x D U L D L Ux D U L D L b D U bC D U L D L U g D U L D L b D U b+--++-------+-----⎧=-+-⎪⎨⎪=--⎩⇒=--+--+-∴=--=--+-分析收敛性：()()()()()1111L D L D L DD LI D D L----=--+-=-+-⎡⎤⎣⎦()()()1111D D L I D L D D D L LD ----⎡⎤=---=-⎣⎦()()111C D U D D L LD U---∴=--()()()()()11112D L D D U L D U D L I DL D U L D U D L U-------=----=--=()()111I C I D L D D U LD U ---⎡⎤∴-=---⎣⎦()()()()111I D L D D U D L D D U A ---⎡⎤⎡⎤=------⎣⎦⎣⎦()()()()1111I I D L D D U A D L D D U A ----⎡⎤⎡⎤=-+--=--⎣⎦⎣⎦ 令()()1M D L D DU -⎡⎤=--⎣⎦1C I M A -∴=-因为A 对称正定，所以D 也正定 令 1111222,()D D D W D D U ----==-TM W W ∴=()()11111112222TWCW D D L LD D D L LD -----⎡⎤⎡⎤=--⎢⎥⎢⎥⎣⎦⎣⎦ 令()11122P DD L LD--=-1T W C W P P -∴=所以C 与T PP 相似，其特征值均为非负实数1111()T I WCW W I C W WM AW W AW -------=-==所以 1I WCW --为对称正定矩阵，其特征值()110WCW λ--C ⇒的特征值()C λ满足()01C λ≤，故该迭代法收敛35.解：1112112111122212Ax Bx b x A Bx A b Bx Ax b x A Bx A b ----⎧+==-+⎧⎪⇒⎨⎨+==-+⎪⎩⎩∴J 法的迭代公式为：(1)()111111(1)()1222110000k k k k J x x A b A B A Bx x A b A B C A B +---+---⎛⎫⎛⎫⎛⎫⎛⎫-∴=+ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭⎝⎭⎛⎫-∴= ⎪-⎝⎭若λ为矩阵1A B --的特征值，对应的特征向量为11111,0n x R x A Bx x λ-∈≠∴-= 11111111111111111111J J x x x A B x C x x x A B x x x x A B x C x x x A B x λλλλλλ----⎡⎤-⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦⎣⎦-⎡⎤⎡⎤⎡⎤⎡⎤===-⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦⎣⎦∴若n 阶矩阵1A B --有特征值12,,,n λλλ，则2n 阶矩阵J C 有特征值12,,,nλλλ±±±38.（1）解：因为系数矩阵A 对称正定，所以可以运用共轭梯度法（CG ）解此方程组 取()()00,0Tx =,()()()()000r p b Ax 0,1T∴==-=-，()()()()()()00r ,r 12p ,Ap α==()()()0,-T10001x xp2α⎛⎫∴=+= ⎪⎝⎭，()()(),0T10003r rAp2α⎛⎫=-=- ⎪⎝⎭,()()()()()(),11r ,r 94r r β==()()(),-T110039p r p 24β⎛⎫=+= ⎪⎝⎭()()()()()()11111r ,r 12p ,Ap α==,()()()()1,-2T2111x x p α=+=()2x 即为所求方程的精确解。
数值分析上机实践题3-2013
数值分析上机实践题第三次上机题目(二分法)第一组:组长：李龙宇，组员：杜彦霖，胡朋，黄湘云，雷盛华，李伟元 用二分法求方程010423=-+x x , ]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第二组: 组长：王宇彬，组员：马泽川，权涛涛，师楠？？，路世伦，仲晓磊 用二分法求方程04442234=++--x x x x ,]0,2[-∈x 的近似根,要求根精确到 510- ,并求二分次数.第三组: 组长： 薛原 ，组员：谢胜权，杨帆，王正奇，肖特，张锡云 用二分法求方程04442234=++--x x x x ,]2,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第四组: 组长：柴春晓 ，组员： 韩静兰，李金慧，刘从，马超群，孟凯悦 用二分法求方程04442234=++--x x x x ,]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第五组: 组长：龙纯鹏，组员：代喜，白鑫，鲍亚强，周邦安，张佳伟 用二分法求方程02=--x x ,]1,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第六组: 组长：何关瑶 ，组员：纪伟亮，侯佳意，李济言，李振华，马文磊 用二分法求方程06cos 2=-++-x e e x x ,]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第七组: 组长：杨钦 ，组员： 王凌宇，吴凯杰，薛小龙，袁权炜，师俊峰 用二分法求方程0232=-+-x x e x ,]1,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第八组: 组长：汪芳 ，组员：张学利，周幸茹，李雨珏，张飞用二分法求方程05.0cos 2=++x x π,]5.1,5.0[∈x 的近似根, 要求根精确到510- ,并求二分次数.第九组: 组长：刘永鸿 ，组员：黄尚政，李超，郭新磊，何奎奎用二分法求 15 的近似根,要求根精确到 510- ,并求二分次数.第十组: 组长：杨吉望 ，组员：龙力，任金雄，王亮，王文强，谢丁波 用二分法求方程 325 的近似根,要求根精确到 510- ,并求二分次数.第十一组: 组长： 张国强，组员： 赵奇，袁硕，郭凯旋，于沛生，鲍宏雷 用二分法求方程 ,05.0c o s 2=++x x π在]2,0[内的近似根,要求根精确到 510- ,并求二分次数.第十二组: 组长：苏映雪 ，组员： 邓晓庆，钟桂平，崔楚轩，高鹏程 用二分法求方程 0797*******=-+--x x x x 的靠近x=2的近似根,要求根精确到 510- ,并求二分次数.备用题:第一组:用二分法求方程 016=--x x , ]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第二组:用二分法求方程 0t a n =-x x ,]5.4,4[∈x 的近似根,要求根精确到 510- ,并求二分次数.补充知识MATLAB 中自带的求根函数:1. roots :求解多项式P(x)=0的根可以用此语句, 输入多项式P(x)的系数(按降幂排列), 输出为P(x)=0的全部根;例如:要求013178)(39=+-+=x x x x P 的根,可以用以下语句:>> fa =[8,0,0,0,0,0,17,0,-3,1]>> gen= roots(fa)运行后输出全部根.2. fsolve: 求解超越方程f(x)=0的根可以用此语句(也可以解多项式方程,但计算量较大), 输入多项式P(x)的系数(按降幂排列), 输出为P(x)=0的全部根调用格式: X = fsolve(F,X0)其中输入函数F(x)的M文件名和解X的初始值X0，X0可以是矩阵或向量。
数值分析上机实验报告
数值分析上机实验报告摘要：本报告是对数值分析课程上机实验的总结和分析，涵盖了多种算法和数据处理方法，通过对实验结果的分析，探究了数值计算的一般过程和计算的稳定性。
1. 引言数值计算是数学的一个重要分支，广泛应用于物理、金融、工程等领域。
本次实验是对数值分析课程知识的实际应用，通过上机实现算法，探究数值计算的可靠性和误差分析。
2. 实验方法本次实验中，我们实现了多种算法，包括：（1）牛顿迭代法求方程的根；（2）高斯消元法求线性方程组的解；（3）最小二乘法拟合数据点；（4）拉格朗日插值法估计函数值；（5）梯形公式和辛普森公式求积分近似值。
对于每个算法，我们都进行了多组数值和不同参数的实验，并记录了相关数据和误差。
在实验过程中，我们着重考虑了算法的可靠性和计算的稳定性。
3. 实验结果与分析在实验中，我们得到了大量的实验数据和误差分析，通过对数据的展示和分析，我们得到了以下结论：（1）牛顿迭代法求解非线性方程的根能够对算法的初始值和迭代次数进行适当的调整，从而达到更高的稳定性和可靠性。
（2）高斯消元法求解线性方程组的解需要注意到矩阵的奇异性和精度的影响，从而保证计算的准确性。
（3）最小二乘法拟合数据点需要考虑到拟合的函数形式和数据的误差范围，采取适当的数据预处理和拟合函数的选择能够提高计算的准确性。
（4）拉格朗日插值法估计函数值需要考虑到插值点的选择和插值函数的阶数，防止出现龙格现象和插值误差过大的情况。
（5）梯形公式和辛普森公式求积分近似值需要考虑到采样密度和拟合函数的选择，从而保证计算的稳定性和收敛速度。
4. 结论通过本次实验的分析和总结，我们得到了深入的认识和理解数值计算的一般过程和算法的稳定性和可靠性，对于以后的数值计算应用也提供了一定的指导和参考。
数值分析上机题3
数值分析上机题目3实验一1．根据Matlab 语言特点，描述Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法。
2．编写Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法的M 文件。
3．给定2020⨯∈R A 为五对角矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡------------------321412132141412132141412132141412132141213O O O O (1)选取不同的初始向量)0(x 及右端面项向量b ，给定迭代误差要求，分别用编写的Jacobi 迭代法和Gauss-Seidel 迭代法程序求解，观察得到的序列是否收敛？若收敛，通过迭代次数分析计算结果并得出你的结论。
(2)用编写的SOR 迭代法程序， 对于(1)所选取的初始向量)0(x 及右端面项向量b 进行求解，松驰系数ω取1<ω<2的不同值，在5)1()(10-+≤-k k x x 时停止迭代，通过迭代次数分析计算结果并得出你的结论。
实验11、 根据MATLAB 语言特点，描述Jacobi 迭代法，Gauss-Seidel 迭代法和SOR 迭代法。
2、 编写Jacobi 迭代法，Gauss-Seidel 迭代法和SOR 迭代法的M 文件。
Jacobi 迭代法function [x1,k]=GS_2(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0; while norm(x1-x0,1)>10^(-7)&k<100 k=k+1;x0=x1;x1=D\((L+U)*x0+b);endk=kx=x1Gauss-Seidel迭代法function [x1,k]=GS_h(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0; while norm(x1-x0,1)>10^(-7)&k<100 k=k+1;x0=x1;x1=(D-L)\U*x0-D\b;endk=kx=x1SOR迭代法function [x1,k]=SOR_h(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0;w=0.96;while norm(x1-x0,1)>10^(-7)&k<100k=k+1;x0=x1;x1=(D-w*U)\(((1-w)*D+w*L)*x0+w*b);endk=kx=x13、采用Jacobi迭代法，Gauss-Seidel迭代法求解五对角矩阵clear,clcA=diag(3*ones(20,1))+diag((-0.5)*ones(19,1),-1)+diag((-0.5)*ones(19,1 ),1)+diag((-0.25)*ones(18,1),-2)+diag((-0.25)*ones(18,1),2);b=sum(A')';[x1,k1]=Jacob_h(A,b)[x2,k2]=GS_h(A,b)运行结果：两种方法都收敛，k1=27,k2=13。
数值分析上机题Matlab(东南大学)3
0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66 0.68 0.70 0.72
152 139 128 119 110 103 96 90 85 80 76 72 68 65 62 59 56 53 51 49 47 45 43 41 39 38
========================================================================================================================
======================================================================================================================================================================== 习题 3_36 ======================================================================================================================================================================== Omega n x1 x2 x3 x4 x5 x6 x7 x8 x9
-0.71279 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281
数值分析第三次作业解答
数值分析第三次作业解答思考题：1：（a ）对给定的连续函数，构造等距节点上的Lagrange 插值多项式，节点数目越多，得到的插值多项式越接近被逼近的函数。
×;(b) 对给定的连续函数，构造其三次样条函数插值，则节点数目越多，得到的样条函数越接近被逼近的函数。
√(c) 高次的Lagrange 插值多项式很常用。
×(d) 样条函数插值具有比较好的数值稳定性。
√3． 以0.1,0.15,0.2为插值节点，计算()f x = Lagrange 插值多项式 2()P x ， 比较2(0)P 和(0)f ，问定理4.1的结果是否适用本问题？ 解： 构造插值多项式：0122022(0.15)(0.2)()0.050.1(0.1)(0.2)()0.050.05(0.1)(0.15)()0.10.05()()()()(0)0;(0)0.1403x x l x x x l x x x l x P x x x x f P --=⨯--=⨯--=⨯=++==在（0，2）区间，5''''''23()(0.2)118.585458f x x f -=≤=从而，对任意的 '''3()(0,0.2),(0)0.05933!f ξξω∈≤ 不存在'''32()(0,0.2),(0)(0)(0)0.14033!f f P ξξω∈=-=。
演示程序：x=0:0.01:0.2; y=x.^(1/2);plot(x,y,'r')pause,hold onx0=[0.1,0.15 ,0.2]; y0=x0.^(1/2); x=0:0.01:0.2; y1=lagrangen(x0,y0,x); plot(x,y1,'b')5：（a ）求()f x x =在节点123452,0.5,0, 1.5,2x x x x x =-=-=== 的三次样条插值（150M M ==）。
数值研究分析上机题(matlab版)(东南大学)
数值分析上机题(matlab版)(东南大学)————————————————————————————————作者：————————————————————————————————日期：数值分析上机报告第一章一、题目精确值为)11123(21+--N N 。
1) 编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。
2) 编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。
3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？二、通用程序clearN=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0);for a=2:N;Sn1=Sn1+1/(a^2-1); endSn2=single(0);for a=2:N;Sn2=Sn2+1/((N-a+2)^2-1); endfprintf('The value of Sn using different algorithms (N=%d)\n',N); disp('____________________________________________________') fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2);disp('____________________________________________________')三、求解结果Please Input an N (N>1):10^2The value of Sn using differentalgorithms (N=100)____________________________________________________Accurate Calculation0.740049Caculate from large to small0.740049Caculate from small to large0.740050__________________________________四、结果分析有效位数n100 10000 1000000顺序从大到小 6 3 3从小到大 5 6 6可以得出，算法对误差的传播又一定的影响，在计算时选一种好的算法可以使结果更为精确。
中科院研究生院信息工程学院课件数值分析数值分析第三次作业及答案
中科院研究⽣院信息⼯程学院课件数值分析数值分析第三次作业及答案6数值分析第三次作业及答案明当h T 0时，它收敛于原初值问题的准确解 y证：梯形公式为 y n ⼗ yn+—[f(X n ，y n )+f(X n^1，y n 』] h由 f (X,y) = —y= y n+ =y n +2( — y n — yn G=L n = 3 Y y n」訓 /乂⼚％l 2+h <12 +h ⼃丫2. (P202(6))写出⽤四阶经典的龙格⼀库塔⽅法求解下列初值问题的计算公式:y n + =y n + — (k 1 +2k 2 +2k 3 +k 4)=0.2214x n +1.2214y n +0.0214 6飞=3y n ⼼+x n )2)” k 2 =3(y n +0.1k 1〃(1+Xn +0.1) )L s =3(yn +0.1k2”(1+Xn +0.1)k 4 =3(yn +0.2k 3”(1+Xn +0.2) 0 2yn + =y n +〒(k 1 +2k 2 +2k 3 +k 4).1. ( P201 (4))⽤梯形⽅法解初值问题〔爲证明其近似解为y n 偌〕:并证⽤ . f 2-h 1 因 yoT " yn F ⼃.⽤上述梯形公式以步长h 经n 步计算到y n ,故有nh ：=x.X◎ T 茹Jf 2—h \n7 l 2+h ⼃1) ]y =x + y, 0 e x ￡1; ly(0) =1;2)l y \3%+x)，O *1； [y(0)=1.解：令h =0.2k 1 = f (X n , y n )= h k2=f (Xn+；；,yn+-k1)=Xn+- + 2 2 2 h k s = f (X n +；, y n +-k 2)=X n +- +y n +-k 2 =1.11(X n + y n )+0.11 2 2 2 2X n +y nh 1)4h h ??yn +；；k i =1.1(Xn +y n )+0.1 2 2 h . . h .................................. .2 ⼋ 2 J 'k 4 = f(X n +h,y n +hk 3)=X n + h + y n +hk 3 =1.222(X n +y n )+0.2223. (P202(7))证明对任意参数t，下列龙格库塔—公式是⼆阶的:r hy n 卄yn+^g+G)；* K i = f (X n, yj；K2 = f (X n +th, y n +thK i);[K3 =f(Xn+(1—t)h,yn+(1—t)hK i).证：由⼀元函数的泰勒展开有2 '''"y(X nG =y(X n) +hy'(X n)⼸[f x(X n,y(X n)) +f y(X n,y(X n))f(X n,y(X n))]中严h'2 3!⼜由⼆元函数的泰勒展开有y n41 =y n +；2(⼼+K3)=y n +；2[(f(X n,y n) + ￡%区『)也+f y(X n, y n)thf (X n, Y n)⼗。
2014级硕士研究生数值分析上机实习报告(答案)
2014级硕士研究生数值分析上机实习 (第一次)姓名： 学号： 学院：实习题目：分别用二分法和Newton 迭代法求方程02010223=-++x x x 的根. 实习目的：掌握两种解法，体会两种解法的收敛速度.实习要求：用C 程序语言编程上机进行计算，精确到8位有效数字. 报告内容：1. 确定实根的个数以及所在区间.2. 将最后两次计算结果填入下表（保留8位数字）：3. 实习过程中遇到哪些问题？如何解决？有何心得体会？4. 两种解法的计算程序（此页写不下时可以加页）：【二分法】#include <stdio.h>#include <math.h>float getvalue(float x){return x*x*x+2*x*x+10*x-20;}void main(){float a=0,b=2,c;c=(a+b)/2;while(fabs(getvalue(c))>0.00001 && fabs(a-b)>0.00001){if(getvalue(c)*getvalue(b)<0) a=c;if(getvalue(a)*getvalue(c)<0) b=c;c=(a+b)/2;}printf("%0.7f\n",c);}【牛顿迭代法】#include "stdio.h"#include "math.h"main(){float x,f,f1;x=8; // x的初值可为任意值do{f=x*x*x+2*x*x+10*x-20; f1=3*x*x+4*x+10; x=x-f/f1;}while(fabs(f)>0.000001);printf("x=%f,f=%f\n",x,f);}2014级硕士研究生数值分析上机实习(第二次)姓名：学号：学院：实习题目：计算8阶三对角矩阵)tridiag235A的行列式.,=.0(235.0,.1274实习目的：掌握计算行列式的方法.实习要求：首先选择一种算法，然后用C程序语言编程上机进行计算.报告内容：1. 简单描述所采用的算法：2. 计算结果：=A3. 实习过程中遇到哪些问题？如何解决？有何心得体会？4. 写出C语言计算程序（此页写不下时可以加页）：#include<stdio.h>#include<math.h>int main(){int n,i,j,k,m,l,SwarpNum=0;double a[10][10],b,temp,result=1;printf("输入行列式阶数:");scanf("%d",&n);printf("输入各值：\n");for(i=0;i<n;i++){for(j=0;j<n;j++)scanf("%lf",&a[i][j]);}for(k=0;k<n-1;k++){if(a[k][k]==0){for(m=n-1;m>k;m--){if(a[m][k]!=0){for(l=0;l<n;l++){temp=a[k][l];a[k][l]=a[m][l];a[m][l]=temp;}SwarpNum++;break;}}}for(i=k+1;i<n;i++){b=-a[i][k]/a[k][k];for(j=k;j<n;j++)a[i][j]=a[k][j]*b+a[i][j];}printf("\n");}for(i=0;i<n;i++)for(j=0;j<n;j++){if(i==j)result*=a[i][j];}result=pow(-1,SwarpNum)*result; printf("result=%f\n\n",result);return 0;}2014级硕士研究生数值分析上机实习 (第三次)姓名： 学号： 学院：实习题目：分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解线性方程组⎪⎩⎪⎨⎧=++=++=++9.14.35.16.84.22.78.17.27.64.38.91.2z y x z y x z y x 实习目的：感受两种迭代法的收敛速度.实习要求：首先构造收敛的Jacobi 迭代法和Gauss-Seidel 迭代法，然后用C 程序语言编程上机进行求解，初始值均取为0，精确到4位小数.报告内容：1. 写出收敛的Jacobi 迭代法和Gauss-Seidel 迭代法：2. 将最后一次迭代次数k与相应的迭代结果填入下表3. 实习过程中遇到哪些问题？如何解决？有何心得体会？4. C语言计算程序：Jacobi迭代法#include <stdio.h>#include <math.h>int function(float y[3],float x[3]); /*判断是否收敛*/float x[3]={0,0,0},z; /*定义初始向量x*/ int i,j,k,n=3;main(){floata[3][3]={{ 8.6,1.5,3.4},{2.1,9.8,3.4},{2.7,1.8,7.2}}, b[3]={1.9,6.7,2.4 };float y[3],sum;int flag;for (k=0;k<100;k++) /*迭代的次数*/{for(i=0;i<n;i++){sum=0;for(j=0;j<n;j++){if(j!=i)sum=sum+a[i][j]*x[j];}y[i]=(b[i]-sum)/a[i][i]; /*算出该迭代时的y[i]*/}for(i=0;i<n;i++){printf("x%d=%-10.6f",i+1,y[i]);}printf("\n"); flag=function(y,x); /*调用函数function*/if(flag==1) /*结束循环*/break; } } int function(float y[3],float x[3]) /*判断是否收*/ { int flag=0;/*标志主函数中的循环是否要结束*/z=fabs(y[0]-x[0]); for(i=0;i<n;i++) if(z<fabs(y[i]-x[i])) z=fabs(y[i]-x[i]); if(z<10e-6) { flag=1; printf("die dai de ci shu shi k=%d\n",k+1); /*输出得到最后结果迭代的次数*/printf("zui hou de jie guo shi:\n"); for(i=0;i<n;i++) printf("x%d=%-10.6f",i+1,y[i]);/*输出方程组的解*/printf("\n"); } else for(i=0;i<n;i++) /*将y[i]的值赋给x[i]进行下一步的迭代*/ x[i]=y[i]; return (flag); }} Gauss-Seidel 迭代法#include<stdio.h> #include<math.h> #include<iostream.h> #define N 3 double Compare(double a[N],double b[N]) { double c=0; int i;for(i=0;i<=N-1;i++) c+=fabs(a[i]-b[i]); return c; } void Gauss_seidel(double A[N][N],doublex[N],double b[N],double precesion) { int i,j,k; double x2[N],x3[N],sum; for(i=0;i<=N-1;i++) { x2[i]=x[i]; x3[i]=x[i]; } k=1; //k 为迭代次数 while(1) { for(i=0;i<=N-1;i++) {sum=0; for(j=0;j<=N-1;j++) {if(j!=i) sum+=A[i][j]*x2[j]; } x[i]=(b[i]-sum)/A[i][i]; x2[i]=x[i];} //输出每一次迭代的结果 printf("第%d 次迭代:\n",k); printf("x3= "); for(i=0;i<=N-1;i++)printf("%lf ",x3[i]); printf("\n");printf("x= ");for(i=0;i<=N-1;i++)printf("%lf ",x[i]); printf("\n"); //判断是否达到度迭代精 if(Compare(x3,x)<=precesion) { printf("达到迭代精度的方程组的解为:\n"); printf("x= ");for(i=0;i<=N-1;i++) printf("%lf ",x[i]); printf("\n"); break;}else{for(i=0;i<=N-1;i++) x3[i]=x[i];k++;continue; }}}void main(){doubleA[N][N]={{ 8.6,1.5,3.4},{2.1,9.8,3.4},{2.7,1.8,7.2} },x[N]={0},b[N]={ 1.9,6.7,2.4 };Gauss_seidel(A,x,b,1e-10);}2014级硕士研究生数值分析上机实习 (第四次)姓名： 学号： 学院：实习题目：分别用复化梯形公式和复化Simpson 公式计算积分dx x xI x⎰=2e2sin 的近似值n I 和m S .实习目的：体会两种复化求积公式的收敛性与收敛速度.实习要求：用C 程序语言编程上机进行计算，结果要有八位有效数字. 报告内容：1. 写出求n I 和m S 的复化求积公式：2. 写出具有八位有效数字的计算结果： =n ；n I = =m ；m S =3. 实习过程中遇到哪些问题？如何解决？有何心得体会？哈尔滨工业大学（威海）实验报告纸4. C语言计算程序（此页写不下时可以加页）：- 11 -。
清华大学高等数值分析_第三次作业答案
高等数值分析第三章作业参考答案1.考虑线性方程组Ax=b,其中A是对称正定矩阵.用Galerkin原理求解方程K=L=Span(v),这里v是一个固定的向量.e0=x∗−x0,e1=x∗−x1证明(e1,Ae1)=(e0,Ae0)−(r,v)2/(Av,v),(∗)其中r=b−Ax0.v应当取哪个向量在某种意义上是最佳的?证明.令x1=x0+αv,那么r1=r−αAv,e1=e0−αv.由Galerkin原理,有(r1,v)=0,因此α=(r,v)/(Av,v).注意到r1=Ae1,r=Ae,有(Ae1,v)=0.于是(e1,Ae1)=(e0−αv,Ae1)=(e0,Ae1)=(e0,Ae0)−α(e0,Av)=(e0,Ae0)−α(r,v)即(∗)式成立.由(∗)式知当v=e0时, e1 A=0最小,即近似解与精确解的误差在A范数意义下最小,算法一步收敛(但是实际中这个v不能精确找到);在最速下降意义下v=r时最佳.2.求证：考虑线性方程组Ax=b,其中A是对称正定矩阵.取K=L=Span(r,Ar).用Galerkin方法求解,其中r是上一步的残余向量.(a)用r和满足(r,Ap)=0的p向量构成K中的一组基.给出计算p的公式.解.设p=r+αAr,(r,Ap)=0等价于(Ar,p)=0.解得α=−(Ar,r)/(Ar,Ar).(b)写出从x0到x1的计算公式.解.设x1=x0+β1r+β2p,那么r1=r−β1Ar−β2Ap,再由Galerkin原理,有(r1,r)=(r1,p)=0,解得β1=(r,r)/(Ar,r),β2=(r,p)/(Ap,p).(c)该算法收敛吗?解.该算法可描述为：(1)选初始x0∈R n,计算初始残差r0=b−Ax0,ε>0为停机准则；(2)对k=0,1,2,...直到 r k <εαk=−(r k,Ar k) (Ar k,Ar k);p k=r k+αAr k;βk=(r k,r k) (Ar k,r k);γk=(r k,p k) (Ap k,p k);r k+1=r k−βk Ar k−γk Ap k;x k+1=x k+βk r k+γk p k.此算法本质上是由CG迭代一步就重启得到的,所以是收敛的,下面给出证法.设用此算法得到的x k+1=x k+¯p1(A)r k,那么e k+1 A=minp1∈P1e k+p1(A)r k A≤ e k+¯p1(A)r k A= e k−¯p1(A)Ae k A≤max1≤i≤n|˜p(λi)| e k A其中0<λ1≤...≤λn为A的特征值,˜p(t)=1−t¯p1(t)是过(0,1)点的二次多项式.当˜p满足˜p(λ1)=˜p(λn)=−˜p(λ1+λn2)时可使max1≤i≤n|˜p(λi)|达到最小.经计算可得min ˜p max1≤i≤n|˜p(λi)|≤(λ1−λn)2(λ1−λn)2+8λ1λn<1故若令κ=λ1/λn,则e k+1 A≤(κ−1)2κ2+6κ+1e k A,方法收敛.3.考虑方程组D1−F−E−D2x1x2=b1b2,其中D1,D2是m×m的非奇异矩阵.取L1=K1=Span{e1,e2,···,e m},L2= K2=Span{e m+1,e m+2,···,e n}.依次用(L1,K1),(L2,K2)按讲义46和47页公式Az∗=r0r0−Az m⊥LW T AV y m=W T r0x m=x0+V(W T AV)−1W T r0各进行一步计算.写出一个程序不断按这个方法计算下去,并验证算法收敛性.用L i=AK i重复上述各步骤.解.对任意给定x0=x(0)1x(0)2,令r=b−Ax0,V1=[e1,e2,...,e m],V2=[e m,e m+1,...,e n].对L i=K i情形,依次用(L1,K1),(L2,K2)各进行一步计算：(L1,K1)(L2,K2)z(1) 1=V1y1z(2)1=V2y2r0−Az(1)1⊥L1r0−Az(2)1⊥L2(V T1AV1)y1=V T1r0,D1y1=V T1r0(V T2AV2)y2=V T2r0,−D2y2=V T2r0x(1)1=x(1)+V1D−11V T1r0x(2)1=x(2)−V2D−12V T2r0得如下算法：(1)选初始x0∈R n,计算初始残差r0=b−Ax0,ε>0为停机准则；(2)对k=1,2,...直到 r k <ε求解D1y1=r k−1(1:m);求解−D2y2=r k−1(m+1:n);x k=x k−1+V1y1+V2y2;r k=r k−1−AV1y1−AV2y2.收敛性:r k=r k−1−AD−11−D−12rk−1=0−F D−12ED−11rk−1Br k−1算法收敛⇔ρ(B)<1⇔ρ(ED−11F D−12)<1.对L i=AK i情形,依次用(L1,K1),(L2,K2)各进行一步计算:(L1,K1)(L2,K2)z(1) 1=V1y1∈K1z(2)1=V2y2∈K2r0−Az(1)1⊥L1=AK1r0−Az(2)1⊥L2=AK2(V T1A T AV1)y1=V T1A T r0(V T2A T AV2)y2=V T2A T r0(D T1D1+E T E)y1=V T1A T r0(D T2D2+F T F)y2=V T2A T r0x(1) 1=x(1)+(D T1D1+E T E)−1V T1A T r0x(2)1=x(2)+(D T2D2+F T F)−1V T2A T r0得如下算法：(1)选初始x0∈R n,计算初始残差r0=b−Ax0,ε>0为停机准则；(2)对k=1,2,...直到 r k <ε求解(D T1D1+E T E)y1=(A T r k−1)(1:m);求解(D T2D2+F T F)y2=(A T r k−1)(m+1:n);x k=x k−1+V1y1+V2y2;r k=r k−1−AV1y1−AV2y2.收敛性:r k=r k−1−A(D T1D1+E T E)−1(D T2D2+F T F)−1A T rk−1(I−B)r k−1算法收敛⇔ρ(I−B)<1⇔0<λ(B)<2.4.令A=3−2−13−2...............−2−13,b=1...2用Galerkin原理求解Ax=b.取x0=0,V m=W m=(e1,e2,···,e m).对不同的m,观察 b−Ax m 和 x m−x∗ 的变化,其中x∗为方程的精确解.解.对于 b−Ax m 和 x m−x∗ ,都是前n−1步下降趋势微乎其微,到第n步突然收敛。
数值分析上机实验指导书讲解
“数值计算方法”上机实验指导书实验一误差分析实验1.1 (病态问题)实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。
对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。
通过本实验可获得一个初步体会。
数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。
病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等) 。
问题提出：考虑一个高次的代数多项式20p(x) =(x-1)(x-2) (x-20)刑(x-k) (1.1)k=1显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。
现考虑该多项式的一个扰动p(x) x1^0 (1.2)其中；是一个非常小的数。
这相当于是对(1.1 )中x19的系数作一个小的扰动。
我们希望比较(1.1)和(1.2)根的差别，从而分析方程(1.1 )的解对扰动的敏感性。
实验内容：为了实现方便，我们先介绍两个MATLAB函数：“roots”和“ poly”。
u = roots(a)其中若变量a存储n+1维的向量，则该函数的输出u为一个n维的向量。
设a的元素依次为&,a2,…，am，则输出u的各分量是多项式方程a1x n a2x n‘ 亠亠a n x a n1 =0 的全部根；而函数b =p o l y ( v的输出b是一个n+1维向量，它是以n维向量v的各分量为根的多项式的系数。
可见“roots” 和“ poly ”是两个互逆的运算函数。
ess = 0.000000001ve = zeros(1,21);ve(2) =essroots(poly(1:20) ve)上述简单的MATLAB程序便得到（1.2）的全部根，程序中的“esS'即是（1.2）中的: 实验要求：（1）选择充分小的ess,反复进行上述实验，记录结果的变化并分析它们。
《数值分析》上机实验报告
数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
1.1 理论依据：设函数在有限区间[a ，b]上二阶导数存在，且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。
当前后两个的差<=ε时，就认为求出了近似的根。
本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。
1.2 C语言程序原代码：#include<stdio.h>#include<math.h>main(){double x2,f,f1;double x1=1.9; //取初值为1.9do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);}1.3 运行结果：1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论：1.使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。
贵州大学数值分析上机实验答案
数值分析上机实验报告课程名称：数值分析上机实验学院：机械工程学院专业：机械制造姓名：******* 学号：********** 年级：12级任课教师：***老师2012年12月30日一．已知A 与b12.38412 2.115237 -1.061074 1.112336 -0.1135840.7187191.7423823.067813 -2.031743 2.11523719.141823 -3.125432 -1.012345 2.189736 1.563849-0.784165 1.112348 3.123124 -1.061074 -3.125A =43215.5679143.123848 2.0314541.836742-1.056781 0.336993 -1.010103 1.112336 -1.012345 3.12384827.108437 4.101011-3.741856 2.101023 -0.71828 -0.037585 -0.1135842.189736 2.031454 4.10101119.8979180.431637-3.111223 2.121314 1.784317 0.718719 1.563849 1.836742 -3.741856 0.4316379.789365-0.103458 -1.103456 0.238417 1.742382 -0.784165 -1.056781 2.101023-3.111223-0.10345814.7138465 3.123789 -2.213474 3.067813 1.112348 0.336993-0.71828 2.121314-1.103456 3.12378930.719334 4.446782 -2.031743 3.123124 -1.010103-0.0375851.7843170.238417-2.213474 4.44678240.00001[ 2.1874369 33.992318 -25.173417 0.84671695 1.784317 -86.612343 1.1101230 4.719345 -5.6784392]TB ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦=（2）用超松弛法求解Bx=b （取松弛因子ω=1.4，x （0）=0，迭代9次）。
数值分析上机实验报告
一、实验目的通过本次上机实验，掌握数值分析中常用的算法，如二分法、牛顿法、不动点迭代法、弦截法等，并能够运用这些算法解决实际问题。
同时，提高编程能力，加深对数值分析理论知识的理解。
二、实验环境1. 操作系统：Windows 102. 编程语言：MATLAB3. 实验工具：MATLAB数值分析工具箱三、实验内容1. 二分法求方程根二分法是一种常用的求方程根的方法，适用于连续函数。
其基本思想是：从区间[a, b]中选取中点c，判断f(c)的符号，若f(c)与f(a)同号，则新的区间为[a, c]，否则为[c, b]。
重复此过程，直至满足精度要求。
2. 牛顿法求方程根牛顿法是一种迭代法，适用于可导函数。
其基本思想是：利用函数在某点的导数值，求出函数在该点的切线方程，切线与x轴的交点即为方程的近似根。
3. 不动点迭代法求方程根不动点迭代法是一种迭代法，适用于具有不动点的函数。
其基本思想是：从初始值x0开始，不断迭代函数g(x)的值，直至满足精度要求。
4. 弦截法求方程根弦截法是一种线性近似方法，适用于可导函数。
其基本思想是：利用两点间的直线近似代替曲线，求出直线与x轴的交点作为方程的近似根。
四、实验步骤1. 二分法求方程根（1）编写二分法函数：function [root, error] = bisection(a, b, tol)（2）输入初始区间[a, b]和精度要求tol（3）调用函数计算根：[root, error] = bisection(a, b, tol)2. 牛顿法求方程根（1）编写牛顿法函数：function [root, error] = newton(f, df, x0, tol)（2）输入函数f、导数df、初始值x0和精度要求tol（3）调用函数计算根：[root, error] = newton(f, df, x0, tol)3. 不动点迭代法求方程根（1）编写不动点迭代法函数：function [root, error] = fixed_point(g, x0, tol)（2）输入函数g、初始值x0和精度要求tol（3）调用函数计算根：[root, error] = fixed_point(g, x0, tol)4. 弦截法求方程根（1）编写弦截法函数：function [root, error] = secant(f, x0, x1, tol)（2）输入函数f、初始值x0和x1，以及精度要求tol（3）调用函数计算根：[root, error] = secant(f, x0, x1, tol)五、实验结果与分析1. 二分法求方程根以方程f(x) = x^2 - 2 = 0为例，输入初始区间[a, b]为[1, 3]，精度要求tol 为1e-6。
数值分析 第三章上机答案
1、程序：x=-1:0.2:1;f=1./(1+25*x.^2);y=polyval(f,x); %计算原函数每个x 值所对应的函数值 p=polyfit(x,y,3) %对（x ，y ）进行三次拟合并输出三次多项式系数 z=polyval(p,x); %计算拟合后的多项式对应x 点的函数值 plot(x,y,'r',x,z,'b') %画图输出：p =1.8222 1.2090 -0.3619 -0.1140则三次曲线拟合的方程为：1140.03619.02090.18222.1)(233--+=x x x x p2、3次和4次多项式拟合：程序：x=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46];%3次拟合%p3=polyfit(x,y,3)xi=0:0.1:1.0;yi=polyval(p3,xi);subplot(1,2,1);plot(x,y,'*',xi,yi,'r');xlabel('x');ylabel('y');%4次拟合%p4=polyfit(x,y,4)xi=0:0.1:1.0;yi=polyval(p4,xi);subplot(1,2,2);plot(x,y,'*',xi,yi,'r');xlabel('x');ylabel('y');输出：p3 =-6.6221 12.8147 -4.6591 0.9266p4 =2.8853 -12.3348 16.2747 -5.2987 0.9427则3次拟合多项式为：9266.06591.48147.126221.6)(233+-+-=x x x x p 4次拟合多项式为：9427.02987.52747.163348.128853.2)(2344+-+-=x x x x x p另一函数拟合：定义函数：function [C,L]=lagran(x,y)%x 为n 个节点的横坐标所组成的向量，y 为纵坐标所组成的向 %C 为所求的牛顿插值多项式的系数构成w=length(x);n=w-1;L=zeros(w,w);for k=1:n+1V=1;for j=1:n+1if k~=jV=conv(V,poly(x(j)))/(x(k)-x(j));% conv 求积，poly(x)将该多项式的系数赋给向量endendL(k,:)=V;endC=y*L输入命令：x=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46];cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x));yy=polyval(cc,xx);plot(xx,yy,'r');hold on ;plot(x,y,'x');xlabel('x');ylabel('y');x=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];y=[0.94 0.58 0.47 0.52 1.00 2.00 2.46];%y 中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据[C,L]=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);hold onplot(xx,yy,'b',x,y,'.');输出：C =40.6746 -110.2183 110.3671 -57.3264 23.4994 -5.4764 0.9400。
数值分析_第三次上机
4.求f(x)=sin x 在[0,π/2]上的最佳一次逼近多项式。
解：设P 1(x)=a 0+a 1x 是f(x) 的最佳一次逼近多项式，则P 1(x)在[0,π/2]上有三个交错点， 满足0<=x 1<x 2<x 3<=π/2。
由于 [f(x)- P 1(x)]’’=(cos x-a 1)’= -sin x 在[0,π/2]上小于0，定号， 故(cos x-a 1)’在[0,π/2]上单调递减，且仅有一个驻点。
故f(x)- P 1(x)在[0,π/2]上只有一个偏差点x 2，满足[f(x)- P 1(x)]’|x=x2 =cos x 2-a 1=0 (1)。
另外两个偏差点x 1=0 ，x 3=π/2 .于是sin 0-a 0 =sin π/2-a 0-π/2a 1 (2)， sin x 2 –a 0-a 1x 2= -( sin 0-a 0) (3) 由（1）（2）（3）式得：a 1=2/π x 2=arccos 2/π=0.88 a 0=-1.18 所以P 1= -1.18+2/π x 。
6.求f(x)=2x 4+3x 3-x 2+1在[-1,1]上的三次最佳一致逼近多项式。
解：设f(x)的三次最佳一致逼近多项式为P 3(x)，由切比雪夫多项式的极性可得 1/2[f(x)- P 3(x)]=1/8T 4(x)=1/8(8x 4-8x 2+1)所以P 3(x)=f(x)-1/4(8x 4-8x 2+1)= 2x 4+3x 3-x 2+1-2x 4+2x 2-1/4 =3x 3+x 2+3/49.求函数f(x)在指定区间上关于Φ(x)=span{1,x}的最佳平方逼近多项式。
(3)f(x)=cosπx, x ∈[0,1]；(4)f(x)=ln x, x ∈[1,2].解：(3)在[0,1]上，经计算得 d 0= ⎰1)(f dx x =0 ，d 1=⎰1)(x dx x f = -2/π2得到法方程组为a 0+1/2a 1=0 ，1/2a 0+1/3a 1= -2/π2 由上面两式解得 a 0=12/π2 ，a 1= -24/π2所以f(x)=cosπx 在[0,1]上的最佳平方逼近多项式为 S 1*=12/π2 -24/π2 x 。
三次样条插值法《数值分析》上机实验作业
昆明理工大学研究生《数值分析》上机实验作业姓名：学号：专业：一、 课题名称课题七 三次样条插值法1、问题提出 设已知数据如下：求f(x)的三次样条插值函数S(x)。
2、要求（1）满足自然边界条件0)0.1()2.0(=''=''s s ；（2）满足第一类边界条件20271.0)2.0(='s ，55741.1)0.1(='s 。
（3）打印输出用追赶法解出的弯曲向量(0M ,1M ,2M ,3M ,4M )和)1.02.0(i S + （i=0,1,2,3,4,5,6,7,8）的值。
并画出)(x S y =的图形。
二、 班级、姓名、学号三、目的和意义由于航空、造船等工程设计的需要而发展起来所谓样条插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性，而且具有较好的稳定性。
今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越广泛地应用。
其中，尤以三次样条插值函数应用最为广泛，如在高速飞机的机翼形体和船体放样等方面的应用，同时在计算机作图方面更是大有作为。
它能够解决一些既有二阶光滑度，又有二阶连续导数的方程，具有良好的收敛性和稳定性。
1. 通过本次实验进一步了解三次样条插值函数，并通过求解三 弯矩方程组得出曲线函数组；2. 通过MATLAB 编程实现求三次样条插值函数的算法，分别考虑 不同的边界条件，同时用追赶法解出弯曲向量和)1.02.0(i S + （i=0,1,2,3,4,5,6,7,8）的值。
四、计算公式首先我们利用)(x S 的二阶导数值),,2,1,0()(j n j M x S j ==''表达)(x S ,因为在区间]x ,[x 1j j +上)()(j x S x S =是不高于三次的多项式，其二阶导数)(x S ''必是线性函数，所以可表示为：]x ,[x x ,h x x M h x x M (x)S 1j j jj1j j1j j+++∈-+-=''对)(x S ''积分两次并利用1j 1j j j y )(x S ,y )(x S ++==，可定出积分 常数，于是得三次样条表达式。
数值分析上机实验报告
数值分析上机实验报告数值分析上机实验报告一、引言数值分析是一门研究利用计算机进行数值计算的学科。
通过数值分析，我们可以使用数学方法和算法来解决实际问题，例如求解方程、插值和逼近、数值积分等。
本次上机实验旨在通过编程实现数值计算方法，并应用于实际问题中。
二、实验目的本次实验的目的是掌握数值计算方法的基本原理和实现过程，加深对数值分析理论的理解，并通过实际应用提高编程能力。
三、实验内容1. 数值求解方程首先，我们使用二分法和牛顿迭代法分别求解非线性方程的根。
通过编写程序，输入方程的初始值和精度要求，计算得到方程的根，并与理论解进行对比。
2. 数值插值和逼近接下来，我们使用拉格朗日插值和最小二乘法进行数据的插值和逼近。
通过编写程序，输入给定的数据点，计算得到插值多项式和逼近多项式，并绘制出插值曲线和逼近曲线。
3. 数值积分然后，我们使用梯形法和辛普森法进行定积分的数值计算。
通过编写程序，输入被积函数和积分区间，计算得到定积分的近似值，并与解析解进行比较。
四、实验步骤1. 数值求解方程（1）使用二分法求解非线性方程的根。
根据二分法的原理，编写程序实现二分法求解方程的根。
（2）使用牛顿迭代法求解非线性方程的根。
根据牛顿迭代法的原理，编写程序实现牛顿迭代法求解方程的根。
2. 数值插值和逼近（1）使用拉格朗日插值法进行数据的插值。
根据拉格朗日插值法的原理，编写程序实现数据的插值。
（2）使用最小二乘法进行数据的逼近。
根据最小二乘法的原理，编写程序实现数据的逼近。
3. 数值积分（1）使用梯形法进行定积分的数值计算。
根据梯形法的原理，编写程序实现定积分的数值计算。
（2）使用辛普森法进行定积分的数值计算。
根据辛普森法的原理，编写程序实现定积分的数值计算。
五、实验结果与分析1. 数值求解方程通过二分法和牛顿迭代法，我们成功求解了给定非线性方程的根，并与理论解进行了对比。
结果表明，二分法和牛顿迭代法都能够较好地求解非线性方程的根，但在不同的问题中，二者的收敛速度和精度可能会有所差异。
数值分析上机报告
数值分析上机实践报告目录第一章标准题部分 (3)1.1 第一题 (3)1.1.1解题的理论依据、算法进行分析及应用条件 (3)1.1.2 MATLAB计算程序 (4)1..1.3 MATLAB运行结果 (7)1..1.4 结果分析与讨论 (7)1.2 第二题 (7)1.2.1解题的理论依据、算法进行分析及应用条件 (7)1.2.2 MATLAB计算程序 (8)1.2.2 MATLAB运行结果 (9)1.2.3 结果分析与讨论 (9)1.3 第三题 (9)1.3.1解题的理论依据、算法进行分析及应用条件 (9)1.3.2 MATLAB计算程序 (10)1.3.3 MATLAB运行结果 (11)1.3.4 结果分析与讨论 (11)1.4 第四题 (11)1.4.1解题的理论依据、算法进行分析及应用条件 (11)1.4.2 MATLAB计算程序 (12)1.4.3 MATLAB运行结果 (13)1.4.4 结果分析与讨论 (13)第二章自主题部分 (14)钢筋混凝土偏心受压柱配筋设计 (14)2.1 引言 (14)2.2 问题研究 (15)2.3 案例分析 (16)2.4 结论 (18)第一章 标准题部分1.1 第一题用三次样条插值的三弯矩法求f(-0.02)和f(2.56)。
1.1.1 解题的理论依据、算法进行分析及应用条件为使问题一般化，本题的三次样条求解使用了对于任意分化的三弯矩插值法。
若记各节点间距n x x x h j j j ,,2,111⋅⋅⋅=-=--，根据理论分析，三次样条插值函数为：))(6()()(6(6)(6)()(112112111131131-------------+--+-+-=j j j j j j j j j j j j jj j j h x x h M y h x x h M y h x x M h x x M x S其中),(1j j x x x -∈。
而式中的诸j M 可由如下的矩阵方程来确定：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡N N d d d M M M10101111122212μλμλμ这里⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧-⋅⋅⋅=+=-=+=--'==---+='--=------+---+-)1,,2,1(1)](1[6),,(6)(6)(61111111111110000100N j h h h h h h y y h y h d x x x f h y y h y y h h d y h y y h d j j jj j j j j j N N N NN N j j j j j j j j j j j μλμ对于上述矩阵方程，由于其系数矩阵为一三对角阵，故我们用追赶法（Thomas 算法）来求解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
CHINA UNIVERSITY OF PETROLEUM
数值分析上机作业(3)
姓名：张立博
学号： 2011213095
院系：化学工程学院
一. 任务：
用MATLAB 语言编写连续函数最佳平方逼近的算法程序（函数式M 文件）。
并用此
程序进行数值试验，写出计算实习报告。
二. 程序功能要求：
在书中Page355或Page345的程序leastp.m （见附一）的基础上进行修改，使其更加完善。
要求算法程序可以适应不同的具体函数，具有一定的通用性。
所编程序具有以下功能： 1. 用Lengendre 多项式做基，并适合于构造任意次数的最佳平方逼近多项式。
可利用递推关系 0112()1,
()()(21)()(1)()/2,3,.....
n n n P x P x x
P x n x P x n P x n n --===---⎡⎤⎣⎦=
2. 被逼近函数f(x)不用内联函数构造，而改用M 文件建立数学函数。
这样，此程序可通过修
改建立数学函数的M 文件以适用不同的被逼近函数（要学会用函数句柄）。
3. 要考虑一般的情况]1,1[],[)(+-≠∈b a x f 。
因此，程序中要有变量代换的功能。
4. 计算组合系数时，计算函数的积分采用变步长复化梯形求积法（见附三）。
5. 程序中应包括帮助文本和必要的注释语句。
另外，程序中也要有必要的反馈信息。
6. 程序输入：（1）待求的被逼近函数值的数据点0x (可以是一个数值或向量)
（2）区间端点：a,b 。
7. 程序输出：（1）拟合系数：012,,,...,n c c c c
（2）待求的被逼近函数值
00001102200()()()()()n n s x c P x c P x c P x c P x =++++
三：数值试验要求：
1. 试验函数：()co s ,[0,4]f x x x x =∈+；也可自选其它的试验函数。
2. 用所编程序直接进行计算，检测程序的正确性，并理解算法。
3. 分别求二次、三次、。
最佳平方逼近函数）x s (。
4. 分别作出逼近函数）x s (和被逼近函数)(x f 的曲线图进行比较。
（分别用绘图函数plot(0x ,s(0x ))和fplot(‘x cos x ’,[x 1 x 2,y 1,y 2])） 四：计算实习报告要求：
1.简述方法的基本原理，程序功能，使用说明。
2.程序中要加注释。
3.对程序中的主要变量给出说明。
4.附源程序及计算结果。
附：
一、 参考程序
Lengendre 多项式作基的函数最佳平方逼近算法程序LEASTP.m
（此程序只适用于对函数x
xe x f =)( 构造最佳平方逼近多项式）
function [c,s]=leastp(x)
%LEASTP.m:least-square fitting with legendre polynomials p1=1;
p2=inline(‘x’,’x’);
p3=inline(‘(3*x^2-1)/2’,’x’); pp1=1;
pp2=inline(‘x.^2’,’x’);
pp3=inline( ‘(3*x.^2-1)/2.*(3*x.^2-1)/2’,’x’); fp1=inline(‘x.*exp(x)’,’x’); fp2=inline(‘(x.^2).*exp(x)’,’x’);
fp3=inline(‘(x*exp(x)).*(3*x.^2-1)/2’,’x’); c(1)=quad(fp1,-1,1)/2;
c(2)=quad(fp2,-1,1)/quad(pp2,-1,1); c(3)=quad(fp3,-1,1)/quad(pp3,-1,1);
s=c(1)+c(2)*p2(x)+c(3)*p3(x); %end
二、被逼近函数用M 文件建立
function f=fun9(x) f=1./(1+x.^2);
三、变步长复化梯形求积公式的算法
()11111.,(()())
22.0,2
3.(),3.15.*2
6.,
7.,,2
b a h b a T f a f b h H x a H H f x x x h x b T T h H T T I T I h h T T ε-=-=
+==+
=+=+<=+-<==
=4.若，则转若则，输出,停机。
转2.
121
0[()()]
2
11
222
((21))221,2,n n n n n j b a
T f a f b H T T T b a b a
f a j n n n -=-=
+=+=+
--++=∑
二、程序代码
被逼近函数用M 文件建立，如下：
三、数值试验要求：
1、试验函数：f(x)=xe x ；逼近函数值的数据点x 0 =0.5；分别求二次、三次最佳平方逼近，n=2，3，4；
（2）、分别作出逼近函数）x s (和被逼近函数)(x f 的曲线图进行比较。
二次最佳平方逼近图形比较：
-1-0.8-0.6-0.4-0.200.20.40.60.81三次最佳平方逼近图形比较：
-1-0.8-0.6-0.4-0.200.20.40.60.81
四次最佳平方逼近图形比较：
-1-0.8-0.6-0.4-0.200.20.40.60.81被逼近函数用M文件建立，如下：。