中国实验快堆额定工况下冷热钠池数值分析

　第38卷第2期原子能科学技术

Vol.38,No.2　2004年3月Atomic Energy Science and Technology

Mar.2004

中国实验快堆额定工况下冷热钠池数值分析

许义军,陆道纲,杨红义,杨福昌

(中国原子能科学研究院快堆工程部,北京　102413)

摘要:应用计算软件STAR 2CD 对中国实验快堆(CEFR )正常运行工况中的额定工况进行了三维数值分析,使用多孔介质模型对屏蔽柱的影响进行了模拟,给出了冷热钠池的三维温度场和流场,与已有热工设计进行了比较,并着重分析了浮升力在数值模拟计算中的影响,为事故工况下的设备动态分析及相应的设备力学分析提供了数据。研究结果为CEFR 的优化设计及事故分析提供了参考数据和技术支持。关键词:三维数值模拟;冷热钠池;工况分析

中图分类号:TL331 文献标识码:A 文章编号:100026931(2004)022*******

Numerical Analysis of China Experimental F ast R eactor

Cold and H ot Plenums U nder Normal Condition

XU Y i 2jun ,L U Dao 2gang ,YAN G Hong 2yi ,YAN G Fu 2chang

(China Institute of A tomic Energy ,P.O.Box 275295,Beijing 102413,China )

Abstract :The CFD software STAR 2CD is used to simulate the normal operating conditions of cold and hot plenums in China Experimental Fast Reator (CEFR ).Complex 32dimension model is established by using porous medium method.The temperature and velocity distribu 2tions of cold and hot plenums are given.These results are compared with the data of the thermal designs of CEFR which have been completed.Meanwhile ,the influence of the buoy 2ancy forces is analyzed.The results of the calculation are valuable for the CEFR design and accident analysis.

K ey w ords :32dimension numerical analysis ;cold and hot plenums ;regime analysis

收稿日期:2003205219;修回日期:2003209222

作者简介:许义军(1973—

),男,河北井陉人,工程师,硕士,核能科学与工程专业中国实验快堆(CEFR )的额定工况是指反应堆处于100%功率下的运行状态,是CEFR 正常运行工况中一个重要工况。额定工况下反应堆运行状态研究是对反应堆事故工况研究的基础,同时又对反应堆安全分析、堆容器及其堆内构件的力学应力分析与评价具有重要意义。本文利用计算流体力学(CFD )软件STAR 2CD

对CEFR 的堆本体进行三维建模分析。

1　计算程序STAR 2CD

STAR 2CD 是专门用于分析涉及流动和质

量传递以及热交换问题的商业化CFD 应用程序系统,该软件具有很强的网格构造能力、强大而丰富的数学物理模型求解器、异常轻巧和集

中的程序结构,使它在许多工程问题中得到广泛应用。STAR(simulation of turbulent flow in arbitrary regions)有很多湍流模型可以选用,如零方程模型、k2ε系列模型、K2L模型及L ES模型,以适应不同的工程需要。目前版本主要采用SIMPL E算法,非结构化网格。

2　冷热钠池建模分析

CEFR冷热钠池均位于主容器内,主要根据其空间内介质钠的温度不同而命名。冷热钠池是介质钠在加热前和经过堆芯加热后的主要载体,与中间热交换器、一回路泵、堆芯一起构成CEFR一回路主热传输系统。

211　计算建模

冷热钠池是一非常复杂的热力学系统,结构上完全模拟相当困难,对许多重要的堆内设备,如中间热交换器、事故热交换器及其主泵等的模拟,必须结合工程实际和软件功能进行简化分析。另外,堆内还有数量众多的屏蔽柱和很薄的堆内支承隔板,由于其结构与堆的整体尺寸相差很大,这也将带来网格构造上的困难,需在建模中给予特别考虑。

212　模型及其简化

在尽可能考虑反应堆真实模拟前提下,对模型做了如下简化:

1)由于所研究区域的对称性,计算中选取180°作为计算区域;

2)热钠池内的3层水平热屏蔽内的钠流动极为缓慢,隔板又非常薄,仔细刻画板的形状需耗费相当多的网格,所以,采用加权方法等效模拟成固体;

3)热钠池中的约150根屏蔽柱中,正对一回路钠循环泵的柱子排列紧密,可认为钠在其中无流动,这些柱子按固体等效处理,剩余柱子采用多孔介质模型进行模拟;

4)对一次钠净化系统,只考虑其传热影响。

213　多孔介质模型及其适用性验证

因热钠池中的屏蔽柱尺寸小、数量多,实际模拟需大量时间和网格,计算的经济性差,所以,有必要使用多孔介质模型进行合理简化。该模型考虑大量的固体构件对流动和传热产生的影响,假设固体构件均匀分布于控制体内,用分布阻力和分布热源分别表征固体构件对动量交换和能量交换的影响。其压降和速度的关系

为:-K i u i=

9ξi。其中,K i为多孔介质的穿透率,K i=αi|u|+βi,αi、βi(i=1,2,3)分别为3个相互正交方向的系数,每个方向的α和β值可相同,也可不同;u为多孔区域的表面速度; 9p/9ξi为压降。

模型计算就是求得各个方向上的穿透率值。为得到在实际钠池模拟中屏蔽柱各方向的穿透率值,构造1个模拟屏蔽柱真实流动的小模型和相应的多孔介质小模型(图1)。通过赋予真实流动的小模型以不同速度值,得到相应的穿透率,再通过线性关系式得到系数α、β值;然后,将α、β值带入到相应的多孔介质小模型,通过设置相同的整体几何结构、物性、边界条件、压降取点位置进行计算,最后将得到的压降和穿透率值进行比较,以确认α、β值是否正确。

二者的比较结果列于表1。由于z向压降极小,程序中的系数设为无穷大。

从表1可看出:二者之间符合较好,可满足工程计算需要,并将其应用到实际的钠池模拟中。

表1　屏蔽柱实际流动模拟和多孔介质流动模拟数据比较

T able1　Comparison of flow d ata betw een real simulation and porous medium simulation

方向

屏蔽柱实际流动多孔介质流动

u/(m?s-1)Δp/Pa K/(kg?m-3?s-1)αβΔp/Pa K/(kg?m-3?s-1)

相对偏差1)/%

R向011171591150662041313518185112217-2137

θ向012691172202662041313570182205315-6174

注:1)相对偏差=(多孔介质流动模拟值-屏蔽柱实际流动模拟值)/屏蔽柱实际流动模拟值

611原子能科学技术第38卷

图1　屏蔽柱实际流动模拟(a)

和多孔介质流动模拟(b)网格图

Fig.1　Grids of real flow simulation in radial

shielding(a)and porous flow simulation(b)

214　网格生成

计算中,网格的生成约占计算工作量的60%,所以,网格的数量、质量及类型的选取均很重要。网格生成考虑了诸如网格的光顺性、正交性、结点分布特性等重要原则,采用整体拉伸法,使用单独的六面体网格、局部加密和网格耦合技术,较为详细地刻画出-6140、-91445m高度的隔板和R=2130m的隔板。冷热钠池计算的复杂程度不同,在描述冷池时

,使用网格较少,而在堆芯出口区域则使用网格较多。共计33754个网格,其中,流体网格19730个。计算中使用的网格示于图2。215　计算方法

计算中的压力修正方法为PISO算法,该算法是SIMPL E法的改进方法,收敛速度快,但每步计算时间较长。采用该算法的主要目的是为提高收敛速度,并为后续的瞬态计算做准备。本计算采用标准的k2ε模型,其经验系数为:Cμ=0109,Cε1=1144,Cε2=1192,Cε3= 1144,Cε4=-0133,Pr k=1,Prε=11219。

图2　冷热钠池网格划分

Fig.2　Grids of cold and hot plenums

216　计算边界条件

根据计算区域不同,使用以下边界条件:

1)热钠池入口温度和流量边界条件按燃料区分为3区(表2),出口为中间热交换器入口(压力边界条件);

2)冷钠池入口为中间热交换器的出口,单个流量为77136kg/s,温度为354.3℃,出口为一回路钠循环泵的入口(压力边界条件);

3)钠液面的散热热流密度为1kW/m2;

4)1台钠循环泵支撑冷却系统带走热量536kW,1台事故热交换器带走热量5215kW,一次钠净化系统引出管带走热量1217kW,堆容器冷却系统带走热量1782kW;

5)在-81225m以下,堆芯围桶内侧温度按360℃等温分布;在-81225m到组件上端-6185m处,温度从360到440℃线性增加。

表2　堆芯燃料分区的流量和温度

T able2　Flux and temperature of CEFR core sub area 分区号描述流量/(kg?s-1)温度/℃

1第一和第二流量区123103546104

2第三和第四流量区144105542112

3除燃料区以外的其它区33157419137

3　计算结果及其分析

311　热钠池

热钠池的范围主要包括在-91445m以上、半径在21430m以内的部分,同时也包括-614m以上半径在31800m以内的部分。

图3为3个流量分区内的流速分布。从流动分布看,最大流速在堆芯出口处,为11243

711

第2期许义军等:中国实验快堆额定工况下冷热钠池数值分析

m/s 。从堆芯第一、第二分区流出的钠流体,首

先受到中心测量柱的阻挡,然后,沿着中心测量柱向上流动,一直到达液面,接着,通过屏蔽柱和屏蔽支撑筒的入口窗,最终流向中间热交换器的入口处。从堆芯第三分区流出的钠流体,因流速较低而受到第二分区中较高流速的影响,在屏蔽柱以内的热钠池区域发生较大范围的搅混,产生一流动涡,其流态示于图4a 。从中可明显看到流体在多孔介质内的流动情况

。

图3　3个入口流量区域内流动示意图

Fig.3　Flow vector in the 3core inlets

半径为2143m 圆筒上的入口窗的流动分析结果表明:在上部,钠的流向是对着中间热交换器的,而最下排的流动则流向屏蔽柱一侧,相比较而言,

这部分的流量较小,不会对额定工况下堆的传热产生大的影响,在事故工况下,这部分倒流量会对堆芯的冷却产生较好的影响。在热池外围的流动示于图4b 。由于存在多个散热冷源,如堆容器冷却系统、泵支撑冷却系统、事故热交换器的冷却等,

这些区域的流动均沿着冷却壁面向下行进,在接近-6140m 处的流动均非常微弱。

图5a 为热钠池屏蔽柱内的温度分布。可以看到,在不同区域,存在着较大的温度差异。最高温度为堆芯出口温度,第一、二分区的平均温度为540℃,第三区的温度较低,约为430℃。从堆芯流出的钠,经搅混在未经过12个入口窗前的热钠池中的平均温度约528℃,

图4　屏蔽柱内(a )和外(b )热钠池流动示意图

Fig.4　Scheme of flow vector in hot plenum radial shielding inner (a )and outer (b )

温度分布较均匀。位于-614m 隔板以上的热池钠温存在自下而上的温度梯度,其中,中间热交换器的入口温度为522℃,在接近底部处的最低温度仅为480℃,位于泵支撑冷却系统的边上(图5b )。同时还可看到:在同一高度,泵支承冷却系统周围的流体也存在较大温差(图6);另外,热钠池较大的温差也反映在沿高度的变化上,图7示出了θ=73°下热钠池内不同半径处温度沿高度的变化。312　冷钠池

冷钠池的形状较为简单,流动形式相对来说较为单一。额定功率下,冷池的入口温度为354℃,单个中间热交换器的流量为77136kg/s ,其流动形式示于图8a 。其中,冷钠池出口流速为017870m/s ,温度为35913℃。冷池温度分布受边界条件影响较大,其中,最高温度为36312℃,位于-614m 的隔板附近;最

11原子能科学技术第38卷

图5

　屏蔽柱内(a )和外(b )热钠池温度示意图

Fig.5　Scheme of temperature in hot plenum

radial shielding inner (a )and outer (b )

图6　泵支承外钠温周向温度分布

Fig.6　Circumferential temperature distribution of sodium out of the pump support structure 高度,m :◆———-41350;●———-41265;

■——

—-31625;Ε———-31575

低温度为35415℃,位于对称边界附近(图

8b )。同理,位于-91445m 隔板以下的钠温均约为360℃,流动形式更为微弱。313　与已有的热工设计的比较

同已有的热工设计参数的比较列于表3。从表3可看出:现计算值与已有热工设计值之间的最大差值出现在事故热交换器出口温度

图7　热钠池高度方向上的温度分布

Fig.7　Altitudinal temperature distribution

of sodium in hot plenum

半径,mm :+———139414;◆———1527;■———1891;●———2408;Ε———3065

图8　冷池内钠流动(a )和温度(b )分布图

Fig.8　Scheme of sodium flow vector (a )and temperature (b )in cold plenum

上,现计算值高出已有热工设计值25℃,三层水平热屏蔽底层和上层温度分别偏高13℃和

11第2期许义军等:中国实验快堆额定工况下冷热钠池数值分析

表3　现有计算同已有的热工设计的比较

T able3　Comparison of the calculated d ata

and the thermal design d ata

参数

温度/℃

已有的热工设计值现计算值中间热交换器入口温度516522

屏蔽柱内热池搅混温度530528

屏蔽柱外热池搅混温度516520

三层水平热屏蔽底层温度403390

三层水平热屏蔽上层温度471485

事故热交换器入口温度516520

事故热交换器出口温度478503冷池平均温度360354

偏低14℃,其余相差不大。

314　影响计算的各种因素

CEFR额定工况下虽为强迫循环流动,但因计算中采用的密度随温度而变,又因冷热钠池在高度方向上存在较大温差,所以,浮升力的影响较大,不可忽略。计算中比较了有无浮升力情况下的温度差异,位于热钠池屏蔽柱外不同高度处两个平面的温度比较结果示于图9a。可以看出:在180°范围内,温度差异很大,局部点的温差可达40℃左右。所以,计算中必须考虑浮升力的影响。

同时,计算中如不考虑多孔介质模型,则计算结果也有较大差异。从图9b可看到,图示位置处考虑与不考虑多孔介质模型的计算温差约5℃(图9b)。所以,考虑多孔介质模型有助于提高计算精度。

4　结论

计算中使用多孔介质模型对复杂的工程问题的解决起了很好的作用。计算得到了CEFR 额定工况下冷热钠池内的温度和速度分布,并给出了堆内某些关键部件的温度场,为设备的力学应力分析提供了数据,为进一步的瞬态分析提供了初始条件。

计算结果表明:

热钠池中存在较为复杂的

图9　不同高度处有无浮升力(a)

和有无多孔介质(b)的温度比较

Fig.9　Scheme of sodium temperature along height with or without buoyancy force(a)and with

or without porous medium(b)

a:高度3575mm,+———有浮升力,●———无浮升力;

高度4987mm,◆———有浮升力,■———无浮升力

b:高度3825mm,Ε———无多孔介质,◇———有多孔介质;

高度4712mm,▲———无多孔介质,●———有多孔介质

流动现象,如堆芯出口区域流体的强烈搅混, 21410m半径上入口窗中下排口的回流等。中间热交换器入口处流动和温度分布状况较复杂,且热池内的钠温沿高度方向存在很大温度梯度,浮升力在计算中有重要作用,不可忽略。

参考文献:

[1]　陶文铨1计算传热学的近代进展[M]1北京:科

学出版社,19981256～2781

021原子能科学技术第38卷

数值分析实验报告模板

数值分析实验报告模板篇一：数值分析实验报告(一)(完整) 数值分析实验报告 1 2 3 4 5 篇二：数值分析实验报告实验报告一题目：非线性方程求解摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。利用二分法求解给定非线性方程的根，在给定的范围内，假设f(x,y)在[a,b]上连续，f(a)xf(b) 直接影响迭代的次数甚至迭代的收敛与发散。即若x0 偏离所求根较远，Newton法可能发散的结论。并且本实验中还利用利用改进的Newton法求解同样的方程，且将结果与Newton法的结果比较分析。前言：（目的和意义）掌握二分法与Newton法的基本原理和应用。掌握二分法的原理，验证二分法，在选对有根区间的前提下，必是收

敛，但精度不够。熟悉Matlab语言编程，学习编程要点。体会Newton使用时的优点，和局部收敛性，而在初值选取不当时，会发散。数学原理：对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b) Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式xk?1?xk?f(xk) f'(xk) 产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 xk?1?xk?rf(xk) 'f(xk) 其中r为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。程序设计：本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x);

数值分析实验报告176453

实验报告插值法数学实验室数值逼近算法设计级 ____________________________ 号 ____________________________ 名 _____________________________ 实验项目名称实验室所属课程名称实验类型实验日期

实验概述：【实验目的及要求】本次实验的目的是熟练《数值分析》第二章“插值法”的相关内容，掌握三种插多项式插值，三次样条插值，拉格朗日插值，并比较三种插值方法的优劣。本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并去实现。【实验原理】《数值分析》第二章“插值法”的相关内容，包括：牛顿多项式插值，三次样条插值，拉格朗日插值的相应算法和相关性质。【实验环境】（使用的软硬件）软件： MATLAB 2012a 硬件：电脑型号：联想 Lenovo 昭阳E46A 笔记本电脑操作系统： Win dows 8专业版处理器：In tel ( R Core ( TM i3 CPU M 350 @2.27GHz 2.27GHz 实验内容：【实验方案设计】第一步，将书上关于三种插值方法的内容转化成程序语言，用 MATLA B 现; 第二步，分别用牛顿多项式插值，三次样条插值，拉格朗日插值求解不同的问题。【实验过程】（实验步骤、记录、数据、分析）实验的主要步骤是：首先分析问题，根据分析设计 MATLA 程序，利用程序算出问题答案，分析所得答案结果，再得出最后结论。实验一：已知函数在下列各点的值为试用4次牛顿插值多项式 P 4（ x ）及三次样条函数 S （ x ）（自然边界条件）对数据进行插值。用图给出{（ X i , y i ）, X i =0.2+0.08i , i=0 , 1, 11, 10 } , P 4 （ x ）及 S （ x ）。值方法：牛顿在MATLAB 件中

数值分析实验报告1

ve=zeros(1,21); ve(n)=ess; r=roots(a+ve); -6-m s=max(abs(r-rr)) end end 利用符号函数：（思考题一）a=poly(1:20); y=poly2sym(a); rr=solve(y) n

很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于MATLAB来进行问题的分析。学号：06450210 姓名：万轩实验二插值法

数值分析实验报告

实验一误差分析实验（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p Λ 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（）中19x 的系数作一个小的扰动。我们希望比较（）和（）根的差别，从而分析方程（）的解对扰动的敏感性。实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a Λ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a Λ 的全部根；而函数 poly(v)b = 的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve = ))20:1((ve poly roots +

数值计算实验报告

(此文档为word格式，下载后您可任意编辑修改！) 2012级6班###(学号)计算机数值方法实验报告成绩册姓名:宋元台学号：成绩：

数值计算方法与算法实验报告学期： 2014 至 2015 第 1 学期 2014年 12月1日课程名称: 数值计算方法与算法专业:信息与计算科学班级 12级5班实验编号： 1实验项目Neton插值多项式指导教师：孙峪怀姓名：宋元台学号：实验成绩：一、实验目的及要求实验目的：掌握Newton插值多项式的算法，理解Newton插值多项式构造过程中基函数的继承特点，掌握差商表的计算特点。实验要求： 1. 给出Newton插值算法 2. 用C语言实现算法二、实验内容三、实验步骤(该部分不够填写.请填写附页)

1.算法分析：下面用伪码描述Newton插值多项式的算法： Step1 输入插值节点数n，插值点序列{x(i),f(i)},i=1,2,……,n,要计算的插值点x. Step2 形成差商表 for i=0 to n for j=n to i f(j)=((f(j)-f(j-1)(x(j)-x(j-1-i)); Step3 置初始值temp=1,newton=f(0) Step4 for i=1 to n temp=(x-x(i-1))*temp*由temp(k)=(x-x(k-1))*temp(k-1)形成 (x-x(0).....(x-x(i-1)* Newton=newton+temp*f(i); Step5 输出f（x）的近似数值newton(x)=newton. 2.用C语言实现算法的程序代码 #includeMAX_N) { printf("the input n is larger than MAX_N,please redefine the MAX_N.\n"); return 1; } if(n<=0) { printf("please input a number between 1 and %d.\n",MAX_N); return 1; } printf("now input the (x_i,y_i)i=0,...%d\n",n); for(i=0;i<=n;i++) { printf("please input x(%d) y(%d)\n",i,i);

数值分析实验报告

实验一、误差分析一、实验目的 1．通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； 2．通过上机计算，了解误差、绝对误差、误差界、相对误差界的有关概念； 3．通过上机计算，了解舍入误差所引起的数值不稳定性。二．实验原理误差问题是数值分析的基础，又是数值分析中一个困难的课题。在实际计算中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。同时，由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法的好坏会影响到数值结果的精度。三．实验内容对20,,2,1,0 =n ，计算定积分 ?+=10 5dx x x y n n . 算法1：利用递推公式 151--=n n y n y , 20,,2,1 =n , 取 ?≈-=+=1 00182322.05ln 6ln 51dx x y . 算法2：利用递推公式 n n y n y 51511-= - 1,,19,20 =n . 注意到 ???=≤+≤=10 10202010201051515611261dx x dx x x dx x , 取 008730.0)12611051(20120≈+≈y .: 四.实验程序及运行结果程序一： t=log(6)-log(5);

n=1; y(1)=t; for k=2:1:20 y(k)=1/k-5*y(k-1); n=n+1; end y y =0.0884 y =0.0581 y =0.0431 y =0.0346 y =0.0271 y =0.0313 y =-0.0134 y =0.1920 y =-0.8487 y =4.3436 y =-21.6268 y =108.2176 y =-541.0110 y =2.7051e+003 y =-1.3526e+004 y =6.7628e+004 y =-3.3814e+005 y =1.6907e+006 y =-8.4535e+006 y =4.2267e+007 程序2： y=zeros(20,1); n=1; y1=(1/105+1/126)/2;y(20)=y1; for k=20:-1:2 y(k-1)=1/(5*k)-(1/5)*y(k); n=n+1; end 运行结果：y = 0.0884 0.0580 0.0431 0.0343 0.0285 0.0212 0.0188 0.0169

(完整版)哈工大-数值分析上机实验报告

实验报告一题目：非线性方程求解摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。前言：（目的和意义）掌握二分法与Newton法的基本原理和应用。数学原理：对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。

Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式产生逼近解x*的迭代数列{x k}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。程序设计：本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可，下面给出主程序。二分法源程序： clear %%%给定求解区间 b=1.5; a=0;

%%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序：clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f);

数值分析实验报告总结

数值分析实验报告总结随着电子计算机的普及与发展，科学计算已成为现代科学的重要组成部分，因而数值计算方法的内容也愈来愈广泛和丰富。通过本学期的学习，主要掌握了一些数值方法的基本原理、具体算法，并通过编程在计算机上来实现这些算法。算法算法是指由基本算术运算及运算顺序的规定构成的完整的解题步骤。算法可以使用框图、算法语言、数学语言、自然语言来进行描述。具有的特征：正确性、有穷性、适用范围广、运算工作量少、使用资源少、逻辑结构简单、便于实现、计算结果可靠。误差计算机的计算结果通常是近似的，因此算法必有误差，并且应能估计误差。误差是指近似值与真正值之差。绝对误差是指近似值与真正值之差或差的绝对值；相对误差：是指近似值与真正值之比或比的绝对值。误差来源见表第三章泛函分析泛函分析概要泛函分析是研究“函数的函数”、函数空间和它们之间变换的一门较新的数学分支，隶属分析数学。它以各种学科

如果 a 是相容范数，且任何满足为具体背景，在集合的基础上，把客观世界中的研究对象抽范数范数，是具有“长度”概念的函数。在线性代数、泛函分析及相关的数学领域，泛函是一个函数，其为矢量空间内的所有矢量赋予非零的正长度或大小。这里以 Cn 空间为例， Rn 空间类似。最常用的范数就是 P-范数。那么当P 取1, 2 ,s 的时候分别是以下几种最简单的情形: 其中2-范数就是通常意义下的距离。对于这些范数有以下不等式: 1 < n1/2 另外，若p 和q 是赫德尔共轭指标，即 1/p+1/q=1 么有赫德尔不等式: II = ||xH*y| 当p=q=2时就是柯西-许瓦兹不等式般来讲矩阵范数除了正定性，齐次性和三角不等式之矩阵范数通常也称为相容范数。象为元素和空间。女口：距离空间，赋范线性空间，内积空间。 1-范数: 1= x1 + x2 +?+ xn 2-范数: x 2=1/2 8 -范数: 8 =max oo ，那外，还规定其必须满足相容性: 所以

数值分析2016上机实验报告

序言数值分析是计算数学的范畴，有时也称它为计算数学、计算方法、数值方法等，其研究对象是各种数学问题的数值方法的设计、分析及其有关的数学理论和具体实现的一门学科，它是一个数学分支。是科学与工程计算（科学计算）的理论支持。许多科学与工程实际问题（核武器的研制、导弹的发射、气象预报）的解决都离不开科学计算。目前，试验、理论、计算已成为人类进行科学活动的三大方法。数值分析是计算数学的一个主要部分,计算数学是数学科学的一个分支,它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现。现在面向数值分析问题的计算机软件有：C,C++,MATLAB,Python,Fortran等。 MATLAB是matrix laboratory的英文缩写，它是由美国Mathwork公司于1967年推出的适合用于不同规格计算机和各种操纵系统的数学软件包，现已发展成为一种功能强大的计算机语言，特别适合用于科学和工程计算。目前，MATLAB应用非常广泛，主要用于算法开发、数据可视化、数值计算和数据分析等，除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。本实验报告使用了MATLAB软件。对不动点迭代，函数逼近（lagrange插值，三次样条插值，最小二乘拟合），追赶法求解矩阵的解，4RungeKutta方法求解，欧拉法及改进欧拉法等算法做了简单的计算模拟实践。并比较了各种算法的优劣性，得到了对数值分析这们学科良好的理解，对以后的科研数值分析能力有了极大的提高。

目录序言 (1) 问题一非线性方程数值解法 (3) 1.1 计算题目 (3) 1.2 迭代法分析 (3) 1.3计算结果分析及结论 (4) 问题二追赶法解三对角矩阵 (5) 2.1 问题 (5) 2.2 问题分析（追赶法） (6) 2.3 计算结果 (7) 问题三函数拟合 (7) 3.1 计算题目 (7) 3.2 题目分析 (7) 3.3 结果比较 (12) 问题四欧拉法解微分方程 (14) 4.1 计算题目 (14) 4.2.1 方程的准确解 (14) 4.2.2 Euler方法求解 (14) 4.2.3改进欧拉方法 (16) 问题五四阶龙格-库塔计算常微分方程初值问题 (17) 5.1 计算题目 (17) 5.2 四阶龙格-库塔方法分析 (18) 5.3 程序流程图 (18) 5.4 标准四阶Runge-Kutta法Matlab实现 (19) 5.5 计算结果及比较 (20) 问题六舍入误差观察 (22) 6.1 计算题目 (22) 6.2 计算结果 (22) 6.3 结论 (23) 7 总结 (24) 附录

数值分析实验报告

学生实验报告实验课程名称开课实验室学院年级专业班学生姓名学号开课时间至学年学期

if(A(m,k)~=0) if(m~=k) A([k m],:)=A([m k],:); %换行 end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去end end x=zeros(length(b),1); %回代求解 x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end y=x; format short;%设置为默认格式显示,显示5位（2）建立MATLAB界面利用MA TLAB的GUI建立如下界面求解线性方程组：详见程序。五、计算实例、数据、结果、分析下面我们对以上的结果进行测试，求解：

? ? ? ? ? ? ? ? ? ? ? ? - = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - - - - 7 2 5 10 13 9 14 4 4 3 2 1 13 12 4 3 3 10 2 4 3 2 1 x x x x 输入数据后点击和，得到如下结果：更改以上数据进行测试，求解如下方程组： 1 2 3 4 43211 34321 23431 12341 x x x x ?? ???? ?? ???? ?? ???? = ?? ???? - ?? ???? - ???? ?? 得到如下结果：

数值分析实验报告资料

机电工程学院机械工程陈星星 6720150109 《数值分析》课程设计实验报告实验一函数插值方法一、问题提出对于给定的一元函数)(x f y =的n+1个节点值(),0,1,,j j y f x j n ==。试用Lagrange 公式求其插值多项式或分段二次Lagrange 插值多项式。数据如下：（1 求五次Lagrange 多项式5L ()x ，计算(0.596)f ,(0.99)f 的值。（提示：结果为(0.596)0.625732f ≈, (0.99) 1.05423f ≈）实验步骤：第一步：先在matlab 中定义lagran 的M 文件为拉格朗日函数代码为： function[c,l]=lagran(x,y) w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if(k~=j) v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; end

第二步：然后在matlab命令窗口输入： >>>> x=[0.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382]; >>p = lagran(x,y) 回车得到： P = 121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845 由此得出所求拉格朗日多项式为 p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845 第三步：在编辑窗口输入如下命令： >> x=[0.4 0.55 0.65 0.80,0.95 1.05]; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718 *x-15.0845; >> plot(x,y) 命令执行后得到如下图所示图形，然后 >> x=0.596; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718 *x-15.084 y =0.6257 得到f（0.596）=0.6257 同理得到f（0.99）=1.0542

数值分析实验报告1

实验一误差分析实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（1.1）中19x 的系数作一个小的扰动。我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。实验容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 poly(v)b = 的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve =

数值分析实验报告

实验五解线性方程组的直接方法实验5.1 （主元的选取与算法的稳定性）问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。实验要求：（1）取矩阵?? ? ?? ?? ?????????=????????????????=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。思考题一：（Vadermonde 矩阵）设 ?? ??????????????????????=? ? ? ?????????????=∑∑∑∑====n i i n n i i n i i n i i n n n n n n n x x x x b x x x x x x x x x x x x A 0020 10022222121102001111 ，，其中，n k k x k ,,1,0,1.01 =+=，（1）对n=2,5,8，计算A 的条件数；随n 增大，矩阵性态如何变化？（2）对n=5，解方程组Ax=b ；设A 的最后一个元素有扰动10-4，再求解Ax=b （3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。（4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？相关MATLAB 函数提示： zeros(m,n) 生成m 行，n 列的零矩阵 ones(m,n) 生成m 行，n 列的元素全为1的矩阵 eye(n) 生成n 阶单位矩阵 rand(m,n) 生成m 行,n 列(0,1)上均匀分布的随机矩阵 diag(x) 返回由向量x 的元素构成的对角矩阵 tril(A) 提取矩阵A 的下三角部分生成下三角矩阵

数值分析实验报告1

实验一误差分析实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（1.1）中19x 的系数作一个小的扰动。我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 poly(v)b = 的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve =

数值分析实验报告3

实验报告实验项目名称数值积分与数值微分实验室数学实验室所属课程名称数值逼近实验类型算法设计实验日期班级学号姓名成绩

实验概述：【实验目的及要求】本次实验的目的是熟练《数值分析》第四章“数值积分与数值微分”的相关内容，掌握复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式。本次试验要求编写复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的程序编码，并在MATLAB软件中去实现。【实验原理】《数值分析》第四章“数值积分与数值微分”的相关内容，包括：复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的相应算法和相关性质。【实验环境】（使用的软硬件）软件： MATLAB 2012a 硬件：电脑型号：联想 Lenovo 昭阳E46A笔记本电脑操作系统：Windows 8 专业版处理器：Intel（R）Core（TM）i3 CPU M 350 @2.27GHz 2.27GHz 实验内容：【实验方案设计】第一步，将书上关于复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的内容转化成程序语言，用MATLAB实现；第二步，分别用以上求积公式的程序编码求解不同的问题。【实验过程】（实验步骤、记录、数据、分析）实验的主要步骤是：首先分析问题，根据分析设计MATLAB程序，利用程序算出问题答案，分析所得答案结果，再得出最后结论。实验：用不同数值方法计算积分 (1) 取不同的步长h.分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？ (2) 用龙贝格求积计算完成问题（1）。 (3)用勒让德多项式确定零点，再代入计算高斯公式，使其精度达到10-4 （1）在MATLAB的Editor中建立一个M-文件，输入程序代码，实现复合梯形求积公式的程序代码如下：

数值计算方法实验报告

差值法实验日志实验题目：插值法实验目的： 1．掌握拉格朗日插值、牛顿插值、分段低次插值和样条插值的方法。 2．对四种插值结果进行初步分析。实验要求：（1）写出算法设计思想；（2）程序清单；（3）运行的结果；（4）所得图形；（5）四种插值的比较；（6）对运行情况所作的分析以及本次调试程序所取的经验。如果程序未通过，应分析其原因。实验主要步骤： 1．已知函数) f满足： (x x0.0 0.1 0.195 0.3 0.401 0.5 f(0.39894 0.39695 0.39142 0.38138 0.36812 x ） 0.35206 （1）用分段线性插值；打开MATLAB，按以下程序输入： x0=-5:5; y0=1./(1+x0.^2); x=-5:0.1:5; y=1./(1+x.^2); y1=lagr(x0,y0,x); y2=interp1(x0,y0,x); y3=spline(x0,y0,x);

for k=1:11 xx(k)=x(46+5*k); yy(k)=y(46+5*k); yy1(k)=y1(46+5*k); yy2(k)=y2(46+5*k); yy3(k)=y3(46+5*k); end [xx;yy;yy2;yy3]' z=0*x; plot(x,z,x,y,'k--',x,y2,'r') plot(x,z,x,y,'k--',x,y1,'r') pause plot(x,z,x,y,'k--',x,y3,'r') 回车得以下图形:

（2）拉格朗日插值。创建M 文件，建立lagr 函数： function y=lagr1(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end 新建一个M 文件，输入： x0=[0.0 0.1 0.195 0.3 0.401 0.5]; y0=[0.39894 0.39695 0.39142 0.38138 0.36812 0.35206]; x=0.0:0.01:0.5; y1=lagr1(x0,y0,x); 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

数值分析实验报告2

一、实验名称复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式及自适应辛普森积分。二、实验目的及要求 1. 掌握复合梯形求积计算积分、复合辛普森求积计算积分、龙贝格求积计算积分和自适应辛普森积分的基本思路和步骤. 2. 培养Matlab 编程与上机调试能力. 三、实验环境计算机，MATLAB 软件四、实验内容 1.用不同数值方法计算积分9 4 ln 1 0-=? xdx x 。（1）取不同的步长h 。分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确指比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善。（2）用龙贝格求积计算完成问题（1）。（3）用自适应辛普森积分，使其精度达到10-4。五、算法描述及实验步骤 1.复合梯形公式将区间[a,b]划分为n 等份，分点x k =a+ah,h=(b-a)/h,k=0,1,...,n ，在每个子区间[x k ,x k +1](k=0,1,...,n-1)上采用梯形公式（），得 )]()([2 )(b f a f a b dx x f b a +-≈ ? （） )]()(2)([2)]()([21 1 110b f x f b f h x f x f h T n k k k n k k n ++=+=∑∑-=+-= （） ),(),(12 )(' '2b a f h a b f R n ∈-- =ηη （）其中Tn 称为复合梯形公式，Rn 为复合梯形公式的余项。 2.复合辛普森求积公式将区间[a,b]划分为n 等份，在每个子区间[x k ,x k +1](k=0,1,...,n-1)上采用辛普森公式（），得 )]()2 (4)([6b f b a f a f a b S +++-= （）

数值分析实验报告97389

数值分析实验报告（第二章）实验题目：分别用二分法、牛顿迭代法、割线法、史蒂芬森迭代法求方程 f(x)=(x2+1)(x?1)5=0 的根x=1，观察不同初始值下的收敛性，并给出结论。问题分析：题目有以下几点要求： 1.不同的迭代法计算根，并比较收敛性。 2.选定不同的初始值，比较收敛性。实验原理：各个迭代法简述二分法：取有根区间[a,b]的重点x0，确定新的有根区间[a1,b1]的区间长度仅为[a,b]区间长度的一版。对压缩了的有根区间[a1,b1]重复以上过程，又得到新的有根区间[a2,b2]，其区间长度为[a1,b1]的一半，如此反复，……，可得一系列有根区间，区间收敛到一个点即为根。牛顿迭代法：不动点迭代法的一种特例，具有局部二次收敛的特性。迭代格式为 x n+1=x n?f(x n) f′(x n) ,n=0,1,2,… 割线法：是牛顿法的改进，具有超线性收敛的特性，收敛阶为1.618. 迭代格式为 x n+1=x n? f(x n) (n)(n?1) (x n?x n?1),n=1,2,… 史蒂芬森迭代法：采用不动点迭代进行预估校正。至少是平方收敛的。迭代格式为 y n=φ(x n) z n=φ(y n) x n+1=x n? (y n?x n)2 z n?2y n+x n 这里φ(x)可采用牛顿迭代法的迭代函数。实验内容：

1.写出该问题的f(x)函数代码如下： function py= f(x) syms k; y=(k^2+1)*(k-1)^5; yy=diff(y,k); py(1)=subs(y,k,x); py(2)=subs(yy,k,x); end 2.分别写出各个迭代法的迭代函数代码如下：二分法： function y=dichotomie(a,b,e) i=2; m(1)=a; while abs(a-b)>e t=(a+b)/2; s1=f(a); s2=f(b); s3=f(t); if s1(1)*s3(1)<=0 b=t; else a=t; end m(i)=t; i=i+1; end y=[t,i+1,m]; end 牛顿迭代法： function y=NewtonIterative(x,e) i=2; en=2*e;m(1)=x; while abs(en)>=e s=f(x); t=x-s(1)/s(2); en=t-x; x=t; m(i)=t; i=i+1; end y=[x,i+1,m]; end 牛顿割线法： function y=Secant(x1,x2,e) i=3; m(1)=x1,m(2)=x2; while abs(x2-x1)>=e s1=f(x1); s2=f(x2); t=x2-(x2-x1)*s2(1)/(s2(1)-s1( 1)); x1=x2; x2=t; m(i)=t; i=i+1; end

数值分析实验报告

数值分析实验报告姓名：张献鹏学号：8 专业：冶金工程班级：重冶二班

目录 1 拉格朗日插值 (1) 问题背景 (1) 数学模型 (1) 计算方法 (1) 数值分析 (2) 2 复化辛普森求积公式 (2) 问题背景 (2) 数学模型 (3) 计算方法 (3) 数值分析 (5) 3 矩阵的LU分解 (5) 问题背景 (5) 数学模型 (6) 理论基础 (6) 实例 (6) 计算方法 (6) 小组元的误差 (8) 4 二分法求方程的根 (9) 问题背景 (9) 数学模型 (9) 计算方法 (9) 二分法的收敛性 (10) 5 雅可比迭代求解方程组 (11) 问题背景 (11) 数学模型 (11) 理论基础 (11) 实例 (11)

计算方法 (12) 收敛性分析 (13) 6 Romberg求积法 (13) 问题背景 (13) 数学模型： (14) 理论基础 (14) 实例 (14) 计算方法 (14) 误差分析 (15) 7 幂法 (16) 问题背景 (16) 数学模型 (16) 理论基础 (16) 实例 (17) 计算方法 (17) 误差分析 (18) 8 改进欧拉法 (18) 问题背景 (18) 数学模型 (18) 理论基础 (18) 实例 (18) 数学模型 (19) 误差分析 (20)

1 拉格朗日插值问题背景对于函数 211 )(x x f += ,55≤≤-x 求拉格朗日插值。10=n ，把)(x f 和插值多项式的曲线画在同一张图上进行比较，观察数值积分中的Lagrange 插值。数学模型取等距差值节点x x =-5+10x /n ，x =0，1，…..，n ，构造n 次lagrange 插值多项式： x x =∑ 11 +x x 2 x x =0 x x +1(x ) (x ? x x 2)x x +1′ (x x ) 当n =10时，十次插值多项式L 10(x )以及函数f (x )的图像可以由Matlab 画出。计算方法： function f= f( x ) f=1./(1+x.^2); end function y=Lagrange(x0,y0,x); n=length(x0); m=length(x); for i=1:m z=x(i); s=; for k=1:n p=; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; End

数值分析实验报告

(此文档为word 格式，下载后您可任意编辑修改！) 姓名：袁义平学号：班级：信息与计算科学二班实验一误差分析实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动其中是一个非常小的数。这相当于是对（1.1）中的系数作一个小的扰动。我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）

的解对扰动的敏感性。实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。上述简单的Matlab 程序便得到（1.2）的全部根，程序中的“ess ”即是（1.2）中的。实验要求：（1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数很小，我们自然感觉（1.1）和（1.2）的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？