电磁场的数值计算方法
物理系0702班学生杜星星
指导老师任丽英摘要：数值计算方法是一种研究并解决数学问题数值近似解的方法，广泛运用于电气、军事、经济、生态、医疗、天文、地质等众多领域。本文综述了电磁场数值计算方法的发展历史、分类，详细介绍了三种典型的数值计算方法—有限差分法、有限元法、矩量法, 对每种方法的解题思路、原理、步骤、特点、应用进行了详细阐述, 并就不同方法的区别进行了深入分析, 最后对电磁场数值计算方法的应用前景作了初步探讨。
关键词：电磁场；数值计算；有限差分法；有限元法；矩量法
引言
自从1864年Maxwell 建立了统一的电磁场理论，并得出著名的Maxwell 方程以来，经典的数学分析方法是一百多年来电磁学学科发展中一个极为重要的手段, 围绕电磁分布边值问题的求解国内外专家学者做了大量的工作。在数值计算方法之前, 电磁分布的边值问题的研究方法主要是解析法，但其推导过程相当繁琐和困难,缺乏通用性,可求解的问题非常有限。上个世纪六十年代以来,伴随着电子计算机技术的飞速发展,多种电磁场数值计算方法不断涌现,并得到广泛地应用,相对于解析法而言,数值计算方法受边界形状的约束大为减少,可以解决各种类型的复杂问题。但各种数值计算方法都有一定的局限性,一个复杂的问题往往难以依靠一种单一方法解决，因此如何充分发挥各种方法的优势,取长补短,将多种方法结合起来解决实际问题,即混合法的研究和应用已日益受到人们的关注。本文综述电磁场的数值计算方法，对三种常用的电磁场数值计算方法进行分类和比较。
1电磁场数值计算方法的发展历史
在上世纪四十年代，就有人试探用数值计算的方法来求解具有简单边界的电磁场问题，如采用Ritz法[1]，以多项式在整个求解场域范围内整体逼近二阶偏微分方程在求解域中的解。五十年代，采用差分方程近似二阶偏微分方程，诞生了有限差分数值计算方法，开始是人工计算，后来采用机械式的手摇计算机计算，使简单、直观的有限差分法得到应用和发展，该方法曾在欧、美风行一时。1964年美国加州大学学者Winslow以矢量位为求解变量，用有限差分法在计算机上成
功地解算了二维非线性磁场
[1]，此后有限差分法在工程电磁场计算领域大为发展。
1965年，Winslow首先将有限元法从力学界引入电气工程中，1969年加拿大MeGill大学P. P. Silvester运用有限元法成功地进行了波导的计算[2]；七十年代初，P. P. Silvester和M. V. K. Chari合作将有限元法应用于二维非线性磁场的计算，成功地计算了直流电机、同步电机的恒定磁场。此后有关有限元法探讨的论文越来越多，有限元法运用的范围由静态场到涡流场到辐射场，由线性场到非线性场，由各向同性媒质到各向异性、要考虑磁滞损耗，由工程电磁场到生物电磁场等等。有人认为有限元法是求解工程电磁场的最有效最成功的方法。
有限元法和有限差分法都是求解边值问题的方法，属于微分方程法。对于开区域或要求求解连续分布场量的区域，这类方法就会受到自身的限制。1972年英国卢瑟福实验室的C.W.Trowbridge等人提出了积分方程法的思想，给出了二维、三维场问题的离散形式[2]，由于此种方法只需离散源区，不需考虑边界条件，所以它较好地解决了无界开域场和要求连续计算场量的问题。该方法计算精度高，但计算量很大。该实验室Sinkin等人又在积分方程法基础上提出了边界积分方程法(又称边界元法)，用此解决线性场的计算，计算量大为减小。此后该室的学者们将积分方程与微分方程法结合起来，提出了求解三维静磁场的双标量位法等。
在解决天线辐射场、散射场问题中，矩量法是一个很重要的数值计算方法。1968年R. F. Harrington发表了专著“Field computation by Moment Method”，对散射场、天线辐射场、波导场等方面的问题起了很好的推进作用。
除以上所介绍的方法外，随着电磁场数值分析的不断发展，各种新方法不断涌现，如计算电场的模拟电荷法，最小二乘配点法，求解磁场的模拟电流法，以及计算场的图论模型法，快速Fourier变换法、有限体元法、无网格计算法等等。各种方法互相配合，出现了一些混合方法，如：矩量法—模拟电荷法、模拟电荷法—有限元法、有限元法—边界元法等，有效地解决了一些实际问题。近年来人工神经网络，小波理论[3]等也引入了电磁场的数值计算中，瞬态电磁场计算如时域有限差分法的应用有了长足的发展。总之随着现有的电磁场数值计算方法的不断深入发展、提高和完善，新的方法不断产生。
在电磁场的数值解法不断发展的同时，人们并没有忘记长期以来所运用的解
析方法。解析法计算结果精确，且可以用解析式表达计算结果，受这些特点吸引，
解析法与数值计算方法相结合形成的半解析法应运而生，也成为了一种主流解算
方法，并还在不断发展。
电磁场数值计算方法发展走向成熟的一个重要标志是：成熟的方法越来越多
地应用于工程实际问题中，商业化通用软件包不断出现[4]。一个商业化软件包通
常由下面几部分组成：
???????????网格图形显示
生、节点形成空调剖分、网格自动产模拟化：数、边界条件几何尺寸、材料性能参数据定义：前处理
???????非线性叠代
求解代数方程组成离散方程组系数矩阵形数据处理
?????????算与显示局部场域分布的精细计
显示受力和损耗计算与图形质区含线性媒质和非线性媒场图显示按要求输出计算结果后处理)(
以上三部分中前、后处理占用了软件包语句的90%以上，编程的主要工作量
在此，而数据处理，也就是我们目前正在学习的数值计算方法仅占软件语句的
10%以内，但它却是占用计算机内存量和消耗CPU 时间的主要部分。
2 电磁场数值计算方法的分类
求解电磁问题的最终要求就是获得满足实际条件的Maxwell 方程的解，借助
于计算数学中的数值算法能够得到大多数电磁问题的近似解。数值算法的基本思
想[5]就是把连续变量函数离散化，把微分方程化为差分方程；把积分方程化为有
限和的形式，从而建立起收敛的代数方程组，然后利用计算机技术进行求解。
数值计算方法从求解方程的形式看，主要分为积分方程法和微分方程法两大类。积分方程法主要有矩量法和边界元法，微分方程法主要有有限差分法和有限元法。对两种方程法的比较，如表一所示。
3 几种重要的数值计算方法
3.1 有限差分法
在电磁场数值计算方法中，有限差分法是应用最早的一种方法。有限差分法以其概念清晰，方法简单、直观，有大致固定的处理和计算模式，具有一定的通用性等特点，在电磁场数值分析领域内得到了广泛的应用。
3.1.1 有限差分法的基本原理
有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。
3.1.2 差分与差商
设函数)(x f 的自变量x 有一小增量h x =?，则)(x f 的增量为
)()()(x f h x f x f -+=? （3.1）
)(x f ?为函数)(x f 的一阶差分。当增量h 足够小，差分f ?与微分df 之间的差才
足够小。一阶差分f ?是自变量x 的函数。按式（1）,计算)(x f ?的差分)(2x f ?称
二阶差分，且
)()()(2x f h x f x f ?-+?=? （3.2）
函数)(x f 的一阶导数)('x f 为
()()x x f dx df x f x ??==
'→?lim 0 应用差分，)('x f 可表示为 '()()()()f x f x h f x f x x h
?+-≈=? （3.3） 故)('x f 可表示为差分)(x f ?除以有限小差分x ?的商，称为差商。
同理，函数)(x f 的二阶导数)(''x f 可表示为 222
1()1()()()()()2()()x x x d f df df dx x dx dx
f x h f x f x f x h h h h f x h f x f x h h +?=-?+---??≈-????
+-+-= （3.4） 3.1.3 差分方程的构造
现以二维静态电、磁场泊松方程的第一类边值问题为例，来具体阐明有限差
分法的应用。设具有平行平面场特征的电磁场场域D ,如图1所示,为一由闭合边
界L 所界定的平面域，其定解条件可表述为
()()y x F y u x u y x u ,,22222
=??+??=? ()D y x ∈, （3.5） ()()y x f y x u L ,,= （3.6）
对于所给定的偏微分方程定解问题，应用有限差分法，首先需从网格剖分
着手决定离散点的分布方式。原则上，可以采用任意的网格剖分方式，但这将直
接影响所得差分方程的具体内容，进而影响解题的经济性与计算精度。为简化问
题，通常采用完全有规律的分布方式，这样在每个离散点上就能得出相同形式的
差分方程，有效地提高解题速度，因而经常采用正方形的网格的剖分方式。现即
以这种正方形网格剖分场域D ，也就是说，用分别与y x ,两坐标轴平行的两簇等
距网格线来生成正方形网格，即
ih x x i ==
..)..........2,1,0(±±=i jh y y j == ..)..........2,1,0(±±=j h 为步长，网格线的交点()j i y x O ,称为节
点，这样D 域就离散化为由网格节点标成
的离散点得集合。
对场域D 中节点()j i y x O ,是一典型
节点，它与周围的1，2，3和4点构成一
个对称星型。设这些离散点上待求函数的
近似值记为 ),(0j i u u =，),1(1j i u u +=，)1,(2+=j i u u ，),1(3j i u u -=，)1,(4-=j i u u
则式（6）可近似离散化为
[][]F j i u j i u j i u h j i u j i u j i u h =-+-++-+-+)1,(),(2)1,(1),1(),(2),1(12
2（3.7） 即
F h j i u j i u j i u j i u j i u 2),(4)1,(),1()1,(),1(=--+-++++ （3.8）
若式（6）F =0,则节点O 上函数u 的值等于其四周相邻点函数值的平均。因
为差分方程（7），（8）只出现待求函数u 在点()j i y x O ,及其四个临近点的值，故
称之为五点差分格式[6]，根据差分方程组解出各离散点处的待求函数值。
3.2 有限元法
x
图1 正方形网格划分 1i +1i -j y -j y 1j y +
传统的变分法在20世纪二三十年代为其新型时期，理论上发展很快，各种
变分问题的最后求解都可归结为解尤拉方程的边值问题，然而只有在一些特殊情
况下尤拉方程才能求出精确解，在大多数情况下，尤拉方程的精确解无法求出。
四五十年代，随着计算机的出现，使其在实际应用中逐渐为比较灵活、通用的有
限差分法所替代。但是，有限差分法在理论上没有以变分原理为基础，因而其收
敛性和数值稳定性往往得不到保证。随后发展形成的有限元法正是变分法与有限
差分法相结合的成果，它取长补短地在理论上以变分原理为基础，在具体方法构
造上又利用了有限差分法网格离散化处理的思想[7]。
3.2.1 有限元法的基本原理
有限元法是以变分原理为基础，将要求解的微分方程型数学模型—边值问
题，首先转化为相应的变分问题，即泛函求极
值问题；然后，利用剖分插值将变分问题离散
化为普通多元函数的极值问题，最终归结为一
组多元的代数方程组，求解该方程组，从而获
得边值问题的数值解。
3.2.2 泛函、变分问题简介
在微积分学形成初期，以数学物理问题为
背景，与多元函数的极值问题相对应，已在几
何、力学上提出了若干个求解泛函极值的问题。
如图2中的质点最速降线问题所述，质点A
从定点),(11y x 自由下滑到定点B ),(22y x ，试求
使滑行时间最短的质点下滑轨道)(x y y =。
图示滑行弧段s d 所需时间为 gy x y gy x v s t 2d 12d sec d d 2'+===α 滑行总时间为 x gy y t x y T x y J x x T d 21d )]([)]([2120??'+=== （3.9）
（3.9）式)]([x y J J =不仅取决于积分端点1x 和2x ，而且取决于)(x y y =的选取。
J 取决于)(x y ，所以J 是函数)(x y 的函数，称之为)(x y 的泛函，记作
y 图2 最速降线问题
)]([x y J 。于是所述之最速降线问题，在数学上就归结为研究泛函)]([x y J 的极值
问题，即
21122[()]min ()0()x x J y x x y x y x y ?==???==?
? （3.10） 泛函的极值问题就称为变分问题。对一般问题而言，可导出下列对应于一个自变
量x 、单个函数)(x y 及其导数)(x y '的已知函数
x y y x F y J x x d ),,(][2
1?'= （3.11） 式中F 为x 、y 和y '的已知函数．．．．
。泛函][y J 的自变量不是一般的自变量，而是一个或几个函数所属的函数族)(x y 。在端点1x 和2x 上分别等于给定值的无数．．
个．
函数)(x y 中，仅有一个)(x y 能使定积分][y J 达到极小值，此函)(x y 数称为极值函数。因此，变分问题就在于寻求使泛函达到极值的该极值函数)(x y ，即分析
研究泛函的极值问题。
3.2.3 泛函的变分与尤拉方程
泛函变分问题的经典解法有两种，一种称之为直接解法，另一类是间接解法。
直接解法是直接把泛函的极值问题近似地转化为一般多元函数的极值问题，用有
限维子空间中的函数去逼近无穷维空间中的极值函数，从而近似求得泛函的极
值。间接解法是将变分问题转化为尤拉方程（微分方程）的定解问题，即边值问
题来求解。
以式（3.11）这种最简形式来推导尤拉方程。设函数)(x y 稍有变化，记作
y y δ+，y δ称之为)(x y 的变分，它反映了整个函数的变化量。这样泛函][y J 的
值也应随之变动，相应于变分y δ的泛函增量为
x y y x F y y y y x F y J y y J J x x d )],,(),,([][][21?'-'+'+=-+=?δδδ （3.12）
将（3.12）式由多元函数的泰勒公式展开
x y y F y y y y F y y F y y F y y F J d }])(2)([21]{[2222222?+''
??+''???+??+''??+??=? δδδδδδ +++=J J J 32δδδ
（3.13） 式中作为泛函增量J ?的线性主部为
?''??+??=21d ][x x x y y F y y F J δδδ （3.14）
J δ
称为泛函][y J 的一次变分（简称变分）。而J 2δ、J 3δ……分别是函数变分y δ及其导数y 'δ的二次、三次齐次式……等的积分，依次称为二次变分，三次
变分……令变分问题的解为)(x y y =，且设极值解)(x y y =稍有变动y y δ+，令
)(x y εηδ= （3.15）
式中ε为任意给定的微量实参数，ε值就确定了),(εx y y =函数族中的某一
曲线，进而确定泛函)],([εx y J 之值；而)(x η是定义于区间],[21x x 且满足
()()021==x x ηη齐次边界条件的可微函数。于是泛函()εη+y J =()[]ε,x y J =()
εΦ就成为变量ε的函数，且当0=ε时获极值函数的解。)(εΦ在0=ε时取得极值的
必要条件是 0dx y y F y y F 021x x 0=??????'??+??=Φ'=Φ'?=δδεε)()( （3.16）
利用分部积分，并根据变分与微分顺序可互换原理，(3.16）式可写为
0)]([2
1='
??-???x x ydx y F dx d y F δ （3.17） 由于（3.17）对任意y δ均成立，故有 0)(='??-??y
F dx d y F （3.18） 方程（3.18）就称为泛函（3.11）的极值问题的尤拉方程。
综上所述，有限元法的基本特点是：
（1）离散化过程保持了明显的物理意义。这是因为，变分原理描述了支配
物理现象的物理学中的最小作用原理（如力学中的最小势能原理、静电学中的汤
姆逊定理等）。
（2）优异的解题能力。 与其他数值计算方法相比较，有限元法在适应场域
边界几何形状，以及媒质物理性质变异情况复杂的问题求解上，有突出的优点。
（3）从数学理论意义上讲，有限元法作为应用数学的一个重要分支，很少
有其他方法应用的这样广泛。它使微分方程的解法和理论面目一新，推动了泛函
分析与计算方法的发展。
3.3 矩量法
矩量法，是近年来在天线、微波技术和电磁波散射等方面广泛应用的一种方法。从这些实际问题涉及开域、激励场源分布形态较为复杂等特征出发，矩量法是将待求的积分方程问题转化为一个矩阵方程的问题[7]，借助于计算机，求得其数值解，从而在所得激励源分布的数值解基础上，即可算出辐射场的分布及其波阻抗等特性参数。
3.3.1矩量法的基本原理
先选定基函数对未知函数进行近似展开，代入算子方程，再选取适当的权函数，使在加权平均的意义下方程的余量等于0，由此将连续的算子方程转换为代数方程。原则上，矩量法可用于求解微分方程和积分方程，但用于微分方程时所得到的代数方程组的系数矩阵往往是病态的，故在电磁场问题中主要用于求解积分方程。
3.3.2 加权余量法
设给定边值问题的场方程统一表述为如下的算子方程，即
g f L =)( （3.19）
已知边界条件为
()b s s r u u =1 （3.20） ()b s s r q n u
=??2 （3.21）
其中: L 是线性算子,f 是待求函数, g 是已知的源。
若u 为精确解，则方程(3.19)和边界条件(3.20),(3.21)应该完全满足。但大多数情况下，不能得到u 的精确解，只能通过数值方法进行估计。
构造一个由有限个线性无关函数i N (i =1,2,…,n )所组成的基函数集合{}N ，借以展开待求函数u 的近似解为
{}{}u N u N u T i
n
i i ==∑=1~ （3.22） 将u
~代入式(3.19)中必然存在误差，即 ()g u L u
R -=~~ (3.23)
取一个归属于试探函数的权函数集合{}W ，令
()0~=-?dV g u L W v j
(j =1,2, …,n ) （3.24） 式(3.24)由n 个方程构成的方程组，它等价于人为地强制近似解u
~，使其因不能精确地满足场方程而导致的误差在平均的含义上等于零。按式(3.24)展开，所构成的各种求解积分或微分方程近似解的方法可被统称为加权余量法[8]。因为
按给定权函数j W 展开式的式(3.24)，即意味着余量g u L R -=~对j
W 取矩的一组平衡式，故式(3.24)的构造亦就被称为矩量法。
基于加权余量式(3.24)，进行移项处理，便得
gdV W dV u L W v j
v j ??=~ (j =1,2, …,n ) （3.25） 将式(3.22)代入式(3.25)的左端，有
()dV N L W u u N L W i v
j n
i i n i i i v j ?∑∑?===??? ??11 （3.26） 为了书写方便，令()()i j i v j n i i N L W dV N L W u ,1≡?∑= 和g W gdV W j v
j ,=?
代入(3.26)式，则可写成
()g W dV N L W u j i v j n
i i ,1=?∑= ()n 21j ??=,, （3.27）
这样，即展开成含n 个未知数i u 的n 个方程。若用矩阵形式表示，则有
[]{}{}g u l = （3.28）
综上所述，矩量法的特点是：矩量法将连续方程离散化为代数方程组, 既适用于求解微分方程, 又适用于求解积分方程。它的求解过程简单, 求解步骤统一, 应用起来比较方便，然而需要一定的数学技巧, 如离散化的程度、基函数与权函数的选取, 矩阵求解过程等。另外必须指出的是, 矩量法可以达到所需要的精确度、解析部分简单, 可计算量很大, 即使用高速大容量计算机, 计算任务也很繁琐。
4 电磁场数值计算的应用前景
电磁场数值计算方法近六十年来发展如此之快，除由于从事这方面的科研人员的努力之外，主要是其研究成果迅速被电机、电器、变压器、加速器、微波器
件、计算机磁头等领域采用
[9]，对改善产品性能、降低生产成本，起着越来越大的作用。
众所周知，任何电磁器件，包括国民经济中应用极其广泛的电机与变压器，其能量转换都是通过电磁场来实现的。但传统的电磁产品的设计方式由于客观条件与手段的限制，把场的实际分布参数当作集中参数处理，不可避免地带来相当大的误差。在迫不得已的情况下，只能用模拟、实验等方法处理，其耗费大、周期长，可借鉴的经验不多，而实验的模拟也不都是有条件的，此外，工农业和日常生活中所用电磁装置越来越多，生产和销售竞争剧烈，因此有效的设计方法显然受到重视。
在采用电磁场数值计算以前的任何方法，即使是十分精巧的代数解析方法，也只能适用于特别简单的几何结构以及一些特殊的假定模型，有时模型甚至简化到不能容忍的程度，但除此之外别无它法。因为实际的电磁场问题，特别是电机电磁场，其边界情况十分复杂，加上铁的饱和以及导体中的涡流效应，解析方法是无能为力的。而数值方法可以模拟复杂的形状，可以适应非线性问题，可以进行涡流的分布计算。
电磁场的数值计算与流体力学分析、温度分析、机械应力分析、电磁力分析和生产计划等等有机联系起来，由计算机完成全部设计，构成所谓CAE[10]系统(Computer Aided Engineering)。CAE系统基本可包括设计与制造的全部过程，产品的设计、制造方案、准备零件、草图、计算、成本核算、生产、机床数控和试验、模拟和产品的自动测试都可能列入到CAE系统中。目前CAE系统在西方国家发展很快，估计每年增长40%，尤其在机械行业中，形成了所谓CADMAT即计算机辅助设计、加工与试验。
从长远观点看，电磁场数值计算要发挥更实际的作用，必须与其他多种学科相结合以形成CAE系统，CAE的发展必然给工业结构上带来巨大的变化。
结束语
数值计算是一门计算的艺术，电磁场的数值计算横跨了多个学科，是数学理论、电磁理论和计算机应用能力的完美组合。通过本次论文设计，我对差分、变分、泛函、加权余量法等数学知识有了深刻的认识，掌握了有限差分法、有限元法、矩量法的基本原理，并能进行简单的运算。但若用数值计算方法解决复杂的电磁场问题
，则在很大程度上依赖于数学知识及计算机编程能力，这就需要进一步系统学习相关的理论知识。
参考文献
[1] 文舸一．计算电磁学的进展与展望[J]．电子学报，1995，23(10)：62-69．
[2] 刘圣民．电磁场的数值方法[M]．武汉：华中理工大学出版社，1991：23-45．
[3] 邓东皋，彭立中．小波分析[J]．数学进展，1991， 20(3)： 294-310．
[4] 洪伟．计算电磁学研究进展［J］．东南大学学报（自然科学版），2002，32（3）：335-339．
[5] 倪光正，杨仕友，钱秀英等．工程电磁场数值计算[M]．北京：机械工业出版社，2003：
123-161．
[6] 盛剑霓．工程电磁场数值分析[M]．西安：西安交通大学出版社，1991： 69-75．
[7] 方静，汪文秉．有限元法和矩量法结合分析背腔天线的辐射特性［J］．微波学报，2000，
16(2)：139-143．
[8] 连汉雄．电磁场理论的数学方法[M]．北京：北京理工大学出版社， 1990：10-15．
[9] 张平文，刘法启，张宇．小波函数值的计算[J]．计算数学，1995(3)：173-185．
[10] 楼仁海，符果行，袁敬闳．电磁理论［M］．成都：电子科技大学出版社，1996：73-102．
Electromagnetic numerical method
Department of Physics 0702 Student Du XingXing
Tutor Ren LiYing
Abstract：Numerical calculation, which is widely used in many fields, such as electric, military affairs, economy, ecology, medical treatment, astronomy, geology and so on is a method of getting the numerical approximate solution in the studying and solving of mathematical problems. This paper reviews the development history and classification of the electromagnetic numerical calculation, and explicitly introduces three typical numerical calculation methods --the finite difference method, the finite element method, and the moment method. Besides, it also gives a thorough elaboration of each method from the angles of problem solving train of thought, principle, step, characteristic and application, and makes a deep analysis of the differences between the
distinct methods. And, in the end, this paper also launches a preliminary probe of the application prospect of the electromagnetic numerical calculation.
Key words：Electromagnetic field; Numerical calculation; Finite difference method; Finite element method; Method of moment
致谢
光阴似箭，岁月如梭，不知不觉我即将走完大学生涯的第四个年头，回想这一路走来的日子，父母的疼爱关心，老师的悉心教诲，朋友的支持帮助一直陪伴着我，让我渐渐长大，也慢慢走向成熟。
首先，我要衷心感谢一直以来给予我无私帮助和关爱的老师们，特别是我的导师任丽英老师。谢谢你们这四年以来对我的关心和照顾，从你们身上，我学会了如何学习，如何工作，如何做人。
其次，我还要真诚地感谢物本0702班的各位同学，在这四年当中，你们给予了我很多帮助，在我的学习、工作、生活各个方面，你们给我提出了很多宝贵的建议，我的成长同样离不开你们。
再次，我要感谢我的父母及家人，没有人比你们更爱我，你们对我的关爱让我深深感受到了生活的美好，谢谢你们一直以来给予我的理解、鼓励和支持，你们是我不断取得进步的永恒动力。
（注：素材和资料部分来自网络，供参考。请预览后才下载，期待你的好评与关注！）
数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名：
课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m )
则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。
目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20)
4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28)
数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10）
1.2牛顿法 程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10）
数值分析上机题 第一章 17.（上机题）舍入误差与有效数 设∑=-= N j N j S 2 2 11 ，其精确值为）111-23（21+-N N 。 （1）编制按从大到小的顺序1 -1 ···1-311-21222N S N +++=，计算N S 的通用 程序； （2）编制按从小到大的顺序1 21 ···1)1(111 222-++--+ -=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算210S ,410S ,610S ,并指出有效位数（编制程序时用单精度）； （4）通过本上机题，你明白了什么？ 解： 程序： （1）从大到小的顺序计算1 -1 ···1-311-21222N S N +++= ： function sn1=fromlarge(n) %从大到小计算sn1 format long ; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end （2）从小到大计算1 21 ···1)1(111 2 22 -++--+-= N N S N function sn2=fromsmall(n) %从小到大计算sn2 format long ; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end （3） 总的编程程序为： function p203()
clear all format long; n=input('please enter a number as the n:') sn=1/2*(3/2-1/n-1/(n+1));%精确值为sn fprintf('精确值为%f\n',sn); sn1=fromlarge(n); fprintf('从大到小计算的值为%f\n',sn1); sn2=fromsmall(n); fprintf('从小到大计算的值为%f\n',sn2); function sn1=fromlarge(n) %从大到小计算sn1 format long; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end function sn2=fromsmall(n) %从小到大计算sn2 format long; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end end 运行结果：
第一章 一、题目 设∑ =-= N N j S 2 j 2 1 1，其精确值为)11 123(21+--N N 。 1) 编制按从大到小的顺序1 1 13112122 2-+??+-+-=N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 N=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); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn (N=%d)\n',N); 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('____________________________________________________')
三、结果 从结果可以看出有效位数是6位。 感想：可以得出，算法对误差的传播有一定的影响，在计算时选一种好的算法可以使结果更为精确。从以上的结果可以看到从大到小的顺序导致大数吃小数的现象，容易产生较大的误差，求和运算从小数到大数所得到的结果才比较准确。
数值分析上机题目 1、 分别用不动点迭代与Newton 法求解方程250x x e -+=的正根与负根。 2、 Use each of the following methods to find a solution in [0.1,1] accurate to within 10^-4 for 4326005502002010x x x x -+--= a. Bisection method b. Newton’s method c. Secant method d. Method of False Position e. Muller’s method 3、 应用Newton 法求f （x ）的零点，e=10^-6,这里f （x ）=x-sin （x ）。 再用求重根的两种方法求f （x ）的零点。 4、 应用Newton 法求f （x ）的零点，e=10^-6,f(x)=x-sin(x) 再用Steffensen’s method 加速其收敛。 5、 用Neville’s 迭代差值算法，对于函数2 1 (),11125f x x x = -≤≤+进行lagrange 插值。取不同的等分数n=5,10，将区间[-1,1]n 等分，取等距节点。把f(x)和插值多项式的曲线画在同一张图上进行比较。 6、 画狗的轮廓图 7、 Use Romberg integration to compute the following approximations to ? a 、 Determine R1,1,R2,1,R3,1,R4,1and R5,1,and use these approximations to predict the value of the integral. b 、 Determine R2,2 ,R3,3 ,R4,4 ,and R5,5,and modify your prediction. c 、 Determine R6,1 ,R6,2 ,R6,3 ,R6,4 ,R6,5 and R6,6,and modify your prediction.
昆明理工大学工科研究生《数值分析》上机实验 学院：材料科学与工程学院 专业：材料物理与化学 学号：2011230024 姓名: 郑录 任课教师：胡杰
P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误！未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? （1）用MA TLAB函数“eig”求矩阵全部特征值。 （2）用基本QR算法求全部特征值（可用MA TLAB函数“qr”实现矩阵的QR分解）。解：MA TLAB程序如下： 求矩阵A的特征值： clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果： 求矩阵B的特征值： clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果：
求矩阵错误！未找到引用源。的特征值： clear; 错误！未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误！未找到引用源。) 输出结果： （2）A= 10 7877565861097 5 9 10 第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到： 第二部：[Q1,R1]=qr(A1);A2=R1*Q1
数值计算方法I 上机实验考试题（两题任选一题） １．小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2. A. 建立火箭升空过程的数学模型（微分方程）； B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度． ２．小型火箭初始质量为1200千克，其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k ，火箭升空过程的数学模型为 0)0(,0,01222==≤≤-+?? ? ??-==t dt dx x t t mg T dt dx k dt x d m 其中)(t x 为火箭在时刻t 的高度，m =1200-15t 为火箭在时刻t 的质量，T （=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t 1 (=900/15=60秒)为引擎关闭时刻． 今测得一组数据如下（t ~时间（秒），x ~高度（米），v ~速度（米/秒））： 现有两种估计比例系数k 的方法： 1．用每一个数据(t,x,v )计算一个k 的估计值（共11个），再用它们来估计k 。 2．用这组数据拟合一个k ． 请你分别用这两种方法给出k 的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）．
数值计算方法作业 姓名：李琦 学号：062410124 求 013=--x x 在x=1.5附近的一个根。 一．牛顿下山法： #include
二．加权法 #include
三．单点弦法： #include
《数值计算方法》上机实验报告华北电力大学 实验名称数值il?算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一. 各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程 *对于非线性方程，若已知根的一个近似值，将在处展开成一阶 xxfx ()0, fx ()xkk 泰勒公式 "f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2! 忽略高次项，有 ，fxfxfxxx 0 ()()()，，, kkk 右端是直线方程，用这个直线方程来近似非线性方程。将非线性方程的 **根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkk fx 0 fx 0 0,
解出 fX 0 *k XX,, k' fx 0 k 水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ik fx ()k 八XX, Ikk* fx()k 这就是牛顿迭代公式。 ,2,计算机程序框图:，见, ,3,输入变量、输出变量说明: X输入变量:迭代初值，迭代精度，迭代最大次数，\0 输出变量:当前迭代次数，当前迭代值xkl ,4,具体算例及求解结果: 2/16 华北电力大学实验报吿 开始 读入 l>k /fx()0?,0 fx 0 Oxx，，01* fx ()0 XX,,，?10 kk, ,1，kN, ?xx, 10 输出迭代输出X输出奇异标志1失败标志
,3,输入变量、输出变量说明: 结束 例：导出计算的牛顿迭代公式，并il ?算。（课本P39例2-16） 115cc （0）, 求解结果: 10. 750000 10.723837 10. 723805 10. 723805 2、列主元素消去法求解线性方程组，1,算法原理: 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角 3/16 华北电力大学实验报告方程组求解。 列选主元是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝 aakkkkkk 对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的 两行（包括常ekk 数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结 ,2,计算机程序框图:，见下页, 输入变量:系数矩阵元素，常向量元素baiji 输出变量:解向量元素bbb,,12n
1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2;
数值分析上机题1 设2 21 1N N j S j ==-∑ ，其精确值为1311221N N ??-- ?+?? 。 （1）编制按从大到小的顺序222 111 21311 N S N = +++---，计算N S 的通用程序。 （2）编制按从小到大的顺序22 21111(1)121 N S N N =+++----，计算N S 的通用程序。 （3）按两种顺序分别计算210S ，410S ，610S ，并指出有效位数。（编制程序时用单精度） （4）通过本上机题，你明白了什么？ 程序代码（matlab 编程）： clc clear a=single(1./([2:10^7].^2-1)); S1(1)=single(0); S1(2)=1/(2^2-1); for N=3:10^2 S1(N)=a(1); for i=2:N-1 S1(N)=S1(N)+a(i); end end S2(1)=single(0); S2(2)=1/(2^2-1); for N=3:10^2 S2(N)=a(N-1); for i=linspace(N-2,1,N-2) S2(N)=S2(N)+a(i); end end S1表示按从大到小的顺序的S N S2表示按从小到大的顺序的S N 计算结果
通过本上机题，看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。从大到小的顺序计算得到的结果的有效位数少。计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。
东南大学数值分析上机作业 汇总 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII
数值分析上机报告 院系： 学号： 姓名：
目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) （1）最大δ值文件 (4) （2）验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序： (10) 2、数据输入及计算结果 (12)
作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ，其精确值为?? ? ??---1112321N N . （1）编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ，计算N S 的通用程序； （2）编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么？ 程序： 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果（省略>>）： S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果（省略>>）： S=cddx(80) S= 0.737577 3、两种方法有效位数对比
实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和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);
第二章 1. 题目：运用MATLAB编程实现牛顿迭代 2. 实验操作 1、打开MATLAB程序软件。 2、在MATLAB中编辑如下的M程序。 function [p1,err,k,y]=newton(f,df,p0,delta,max) %f 是要求根的方程(f(x)=0); %df 是f(x)的导数； %p0是所给初值,位于x*附近； %delta是给定允许误差； %max是迭代的最大次数； %p1是newton法求得的方程的近似解； %err是p0的误差估计； %k是迭代次数； p0 for k=1:max p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; k p1 err y=feval('f',p1) if (err
p0 = 1.2000 k =1 p1=1.1030 err=0.0970 y=0.0329 k= 2 p1=1.0524 err=0.0507 y=0.0084 k =3 p1=1.0264 err=0.0260 y=0.0021 k =4 p1=1.0133 err=0.0131 y=5.2963e-004 k =5 p1=1.0066 err=0.0066 y=1.3270e-004 k =6 p1=1.0033 err=0.0033 y=3.3211e-005 k =7 p1=1.0017 err=0.0017 y=8.3074e-006 k =8 p1=1.0008 err=8.3157e-004 y = 2.0774e-006 k =9 p1=1.0004 err=4.1596e-004 y =5.1943e-007 k=10 p1=1.0002 err=2.0802e-004 y= 1.2987e-007 k=11 p1=1.0001 err=1.0402e-004 y =3.2468e-008 k=12 p1=1.0001 err=5.2014e-005 y=8.1170e-009 k=13 p1=1.0000 err=2.6008e-005 y= 2.0293e-009 k=14 p1=1.0000 err=1.3004e-005 y=5.0732e-010 k=15 p1 =1.0000 err=6.5020e-006 y=1.2683e-010 k=16 p1 =1.0000 err=3.2510e-006 y=3.1708e-011 k=17 p1 =1.0000 err=1.6255e-006 y =7.9272e-012 k=18 p1 =1.0000 err =8.1279e-007 y= 1.9820e-012 ans = 1.0000 结果说明：经过18次迭代得到精确解为1，误差为8.1279e-007。
数值计算方法上机实习题 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+ -=-，从I 0=0.1824, 0=0.1823I 出发，计算20I ； （2） 20=0I ，20=10000I , 用n I I n n 51 5111+- =--，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 答：第一个算法可得出 e 0=|I 0?I 0 ?| e n =|I n ?I n ?|=5n |e 0| 易知第一个算法每一步计算都把误差放大了5倍，n 次计算后更是放大了5n 倍，可靠性低。 第二个算法可得出 e n =|I n ?I n ?| e 0=(15 )n |e n | 可以看出第二个算法每一步计算就把误差缩小5倍，n 次后缩小了5n 倍，可靠性高。
2． 求方程0210=-+x e x 的近似根，要求41105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 计算根与步数程序： fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x; f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1); fprintf('root=%6.8f ,n=%d \n',root,n); 计算结果显示： root=0.09057617 ,n=11 （2） 取初值00=x ，并用迭代10 21 x k e x -=+；
（3） 加速迭代的结果; （4） 取初值00 x ，并用牛顿迭代法；
第四章 问题一 一、问题综述 在离地球表面高度为y处的重力加速度如下： 计算高度y=55000m处的重力加速度值。 二、问题分析 以高度y作为自变量，重力加速度的值为因变量。得到以下信息： f(0)=9.8100; f(30000)=9.7487; f(60000)=9.6879; f(90000)=9.6278; f(120000)=9.5682; 本题要求的就是f(55000)的值。 以下将采用课堂中学到的Lagrange插值多项式法、Newton插值多项式法、分段低次插值法和样条插值法求解该问题。 三、问题解决 1. lagrange插值多项式法 对某个多项式函数，已知有给定的k+ 1个取值点： 其中对应着自变量的位置，而对应着函数在这个位置的取值。 假设任意两个不同的x j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：
其中每个为拉格朗日基本多项式（或称插值基函数），其表达式为： 拉格朗日基本多项式的特点是在上取值为1，在其它的点上取值为0。 源程序lagrange.m function [c,f]=lagrange(x,y,a) % 输入：x是自变量的矩阵；y是因变量的矩阵；a是要计算的值的自变量； % 输出：c是插值多项式系数矩阵；f是所求自变量对应的因变量； m=length(x); l=zeros(m,m); % l是权矩阵 f=0; for i=1:m v=1; for j=1:m if i~=j v=conv(v,poly(x(j)))/(x(i)-x(j)); % v是l_i(x)的系数矩阵 end end l(i,:)=v; % l矩阵的每一行都是x从高次到低次的系数矩阵 end c=vpa(y*l,10); % 对应阶次的系数相加，乘以y，显示10位有效数字 for k=1:m f=f+c(k)*a^(m-k); end 输入矩阵 x=[0 30000 60000 90000 120000] y=[9.81 9.7487 9.6879 9.6278 9.5682] a=55000 再运行源函数，可得： c = [ -2.057613169e-23, 4.938271605e-18, -3.703703702e-14, -0.000002046111111, 9.81] f = 9.6979851723251649906109417384537
如有帮助欢迎下载支持 数值分析上机题 姓名：陈作添 学号： 040816 习题 1 20．（上机题）舍入误差与有效数 N 1 1 3 1 1 设 S N ，其精确值为 。 2 2 2 N N 1 j 2 j 1 （1）编制按从大到小的顺序 1 1 1 ，计算 S 的通用程序。 S N 1 32 1 N 2 1 N 2 2 （2）编制按从小到大的顺序 1 1 1 ，计算 S 的通用程序。 S N 1 (N 1)2 1 22 1 N N 2 （3）按两种顺序分别计算 S 102 ， S 104 ， S 106 ，并指出有效位数。 （编制程序时用单精度） （4）通过本上机题，你明白了什么？ 按从大到小的顺序计算 S N 的通用程序为： 按从小到大的顺序计算 S N 的通用程序为： #include
2015.1.9 上机作业题报告 JONMMX 2000
1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果
1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序
数值计算大作业 题目一、非线性方程求根 1.题目 假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。 （1）如果令()N t 表示在t 时刻的人口数目，β 表示固定的人口出生率，则人口数目满足微分方程() ()dN t N t dt β=，此方程的解为0()=t N t N e β； （2）如果允许移民移入且速率为恒定的v ，则微分方程变成() ()dN t N t v dt β=+， 此方程的解为 0()=+ (1) t t v N t N e e βββ -； 假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率β，精确到 410-；且通过这个数值来预测第二年年末的人口数，假设移民速度v 保持不变。 435000 1564000=1000000(1) e e βββ + - 2.数学原理 采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程0)(=x f ，如果) (x f 是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。 设已知方程0)(=x f 有近似根k x (假定0)(≠'x f )，将函数)(x f 在点k x 进行泰勒展开，有 . ))(()()(???+-'+≈k k k x x x f x f x f 于是方程0)(=x f 可近似地表示为 ))(()(=-'+k k x x x f x f 这是个线性方程，记其根为1k x +，则1k x +的计算公式为
数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出： 考虑一个高次的代数多项式 ∏=-= ---=20 1)()20)...(2)(1()(k k x x x x x p （E1-1） 显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简 单的）。现考虑该多项式方程的一个扰动 0)(19 =+x x p ε (E1-2) 其中ε是一个非常小的数。这相当于是对（E1-1）中19 x 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。 实验内容： 为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数 u ＝roots （a ） 其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a 的全部根，而函数 b=poly(v) 的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“Poly ”是两个互逆的运算函数. ve=zeros(1,21); ve(2)=ess; roots(poly(1:20))+ve) 上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（E1-2）中的扰动项改成18 x ε或其他形式，实验中又有怎样的现象出现？