青花漫舞6305467
范文一：数值计算(数值分析)试题及答案
武汉理工大学研究生课程考试标准答案
用纸
课程名称：数值计算（A ） 任课教师 ：
一. 简答题，请简要写出答题过程（每小题5分，共30分）
绝对误差和相对误差分别是多少？
3分）
2分）
01
32. 已知f (x )=x 8+x 5-32，求f ??,3,
01
?,38?f 3，??,3,
,39??.
(5分）
3. 确定求积公式
?
1
f (x ) dx ≈A 0f (0)+A 1f (1)+A 2f '(0)中的待定系数，使其代数
精度尽量高，并指明该求积公式所具有的代数精度。
解：要使其代数精度尽可能的高，只需令f (x ) =1, x ,
x m , 使积分公式对尽可能
大的正整数m 准确成立。由于有三个待定系数，可以满足三个方程，即m =2。 由f (x ) =1数值积分准确成立得：A 0+A 1=1 由f (x ) =x 数值积分准确成立得：A 1+A 2=1/2 由f (x ) =x 2数值积分准确成立得：A 1=1/3
解得A 1=1/3, A 2=1/6, A 0=2/3. (3分）
此时，取f (x ) =x 3积分准确值为1/4, 而数值积分为A 1=1/3≠1/4, 所以该求积公式的最高代数精度为2次。 (2分）
?10-1?
?的谱半径。 0104. 求矩阵A =???
??-202??
λ-1
解 λI -A =
01
0=λ(λ-)(1λ-λ- )3
2
λ-1
矩阵A 的特征值为λ1=0, λ2=1, λ3=3 所以谱半径ρ(A )=max {0,1,3}=3 (5分）
?10099?
5. 设A = ?, 计算A 的条件数cond (A )p , (P =2, ∞).
?9998?
?98-99?A *?-9899?-1
解：A = = ??A =?
-9910099-100A ????
*
矩阵A 的较大特征值为198.00505035，较小的特征值为-0.00505035，则
cond (A ) 2=A 2?A -1c o n (d ) A ∞=
2
=198.00505035/0.00505035=39206(2分）
390 1 6 (3分）
A ∞
-1
=199?199=
∞
H 3(x ) =(1-2
x -x 0x -x 12x -x 1x -x 02
)() y 0+(1-2)() y 1
x 0-x 1x 0-x 1x 1-x 0x 1-x 0
x -x 02x -x 12
'+(x -x 1)('+(x -x 0)() y 0) y 1
x 0-x 1x 1-x 0
(5分）
1
并依条件H (0)=1, H '(0)=, H (1)=2, H '(1)=2. ，得
21
H 3(x ) =(1+2x )(x -1) 2+2(3-2x ) x 2+x (x -1) 2+2(x -1) x 2
2
(5分）
131
=x +x +122
2. 已知f (-1)=2, f (1)=1, f (2)=1，求f (x )的Lagrange 插值多项式。
解：注意到：
x 0=-1, x 1=1, x 2=2; y 0=2, y 1=1, y 2=1l 0=l 1=l 2=
(x -x 1)(x -x 2) (x -1)(x -2)
=
(x 0-x 1)(x 0-x 2) 6(x -x 0)(x -x 2) (x +1)(x -2)
=
(x 1-x 0)(x 1-x 2) -2(x -x 0)(x -x 1) (x +1)(x -1)
=
(x 2-x 0)(x 2-x 1) 3
n
2
j j j =0
∴L (x ) =∑y l (x )=(x -1)(x -2) ?2+(x +1)(x -2) ?1+(x +1)(x -1) ?1
6-23
=
12
x -3x +8)(6
3.3. 给出如下离散数据，试对数据作出线性拟合
解: P (x ) =a 0+a 1x
m m
?
ma 0+(∑x i ) a 1=∑y i
??i =1i =1
(5分） ?m m m
2?(∑x i ) a 0+(∑x i ) a 1=∑x i y i
?i =1i =1?i =1
?4a 0+6a 1=12
?
6a +14a =251?0
a 0=0. 9, a 1=1. 4，P (x ) =0. 9+1. 4x (5分）
?20x 1+2x 2+3x 3=24
T ?(0)
4. 用Jacobi 迭代法求解方程组?x 1+8x 2+x 3=12，取初值x =(0,0,0)，计算
?2x -3x +15x =30
23?1
迭代二次的值；（2分）
问Jacobi 迭代法是否收敛？为什么？（2分）
若收敛，需要迭代多少次，才能保证各分量的误差绝对值小于10-
6？（提示：
问Gauss-Seidel 迭代法是否收敛？为什么？（1分）
解：先将方程组化成便于迭代的形式，以20,8,15分别除以三个方程两边得
136??x =-x -x +0123 ?10205
?
133? -1x =-x -x +B = ， 迭代矩阵?213 8882? 212? x =-x +x +2-23 ?3
155??15
由于||B ||∞=
-
1
100
-
315
3?20??1?-,
?8?0???
1
=q <1, 或者因为原方程组系数矩阵严格对角占优，故jacobi="">
Gauss-Seidel 迭代法收敛。
由||x *-x
(k )
q k ||≤||x (1)-x (0)||得公式 K >1-q
ln
ε(1-||B ||∞)
||x (1)-x (0)||∞
及||x (1)-x (0)||∞
ln ||B ||∞
1
10-6(1-)
1ln(?10-6) ln
ln(3?106) ==≈13.57 可得K >
ln 3ln ln 33
所以迭代14次时，能保证各分量的误差绝对值小于10.
-6
2
??y '=-x +y
5. 用欧拉法解初值问题?在[0,1.5]上的数值解，取h =0.5，计算过程保
??y (0)=2
留5位小数。（要求写出迭代公式，不写公式扣4分）
解：欧拉法的公式为
y (x k )≈y k =y k -1+hf (x k -1, y k -1)=y k -1+h (-x k -1+y k 2-1)，k =1,2,3,4 （4分）
已知x 0=0，y 0=2
y (0.5)≈y 1=2+0.5?(-0+22)=4 y (1)≈y 2=4+0.5?(-0.5+42)=11.75
y (1.5)≈y 3=11.75+0.5?(-1+11.752)=80.28125 （6分）
三. 分析题，请写出主要分析与认证过程（每小题5分，共10分） 1. 设Ax =b ，其中A ∈R n ?n 为非奇异矩阵，证明cond (A T A )=??cond (A )2?? 2
2
证明： (AAT ) =(AT ) T A T =AA T ?AA T 为对称矩阵
T
对于任意给定的非零列向量X, 都有XT A T A X=(AX) T (AX) =b 2>0
所以A A 为正定矩阵, ?AA 也为正定矩阵所以AA T 为对称正定矩阵.
cond(A T A ) 2=A T A
2
T T
(AT A) -1
2
=ρ(AT A) ρ[(AT A) -1]
=
λmax (AT A) λmax [(AT A) -1]
2
又由于A 2=A -1?cond[(A ) 2]2=(A
2
2
=A -1) 2=λmax (AT A) λmax [(AT A) -1]
2
所以cond(A T A ) 2=cond[(A ) 2] （5分）
证明：设x j 是向量X 的分量，则X
2∞
2
n
2∞
2
=?max x i ?≤∑x i ≤n X ?i ?i =1
，
所以由向量范数的概念可知，结论成立. （5分）
证：因为f '(x ) =-
，故牛顿迭代格式为 x 2
1-a x f (x )
=x k -k =x k (2-ax k ), k =0,1, 2, x k +1=x k -
1f '(x ) -2x k
下证明其收敛性。
记第k 步的误差为e k =x *-x k 和构造r k =ae k ，k =0,1, 则有a , r k , x k 三者之间的关系为
r k =ae k =a (x *-x k ) =ax *-ax k =1-ax k ； 而 r k +1=1-ax k +1=1-ax k (2-ax k ) =1-ax k (1+r k ) =(1-ax k ) -ax k r k =r k -ax k r k =(1-ax k ) r k =r k 2, （+）式是一个递推关系，重复使用它，得 r k =r k 2-1=r k 2-1=若0
2
（5分）
，
k =0,1,
（+）
=r , k =0,1,
2k 0
（*）
2
，那么 -1
也即有 |r 0|<1>
r 02=0, 即lim r k =0。 从而有 l i m
k →∞
k →∞
k
又因为r k =ae k ，所以lim e k =0，
k →∞
也就是牛顿法产生的序列{x k }收敛于x *=
1
。（5分） a
范文二：数值分析_数值计算小论文
Runge-Kutta 法的历史发展与应用
摘要 Runge-Kutta 法是极其重要的常微分方程数值解法,本文仅就其起源及发展脉络加以 简要研究。 对 Runge 、 Heun 以及 Kutta 等人的贡献做出适当评述, 指出 Runge-Kutta 方法起 源于 Euler 折线法。同时对 Runge-Kutta 法的应用做简要研究。
关键词 Euler折线法 标准四阶 Runge-Kutta 法 应用
一、发展历史 [1]
1.1 Euler折线法
在微分方程研究之初, 瑞士数学家 L.Euler(1707.4— 1783.9) 做出了开创性 的工作。 他和其他一些数学家在解决力学、 物理学问题的过程中创立了微分方程 这门学科。 在常微分方程方面, Euler 在 1743年发表的论文中, 用代换 kx y e =给 出了任意阶常系数线性微分方程的古典解法,最早引入了“通解”和“特解”的 概念。
1768年, Euler 在其有关月球运行理论的著作中, 创立了广泛用于求初值问 题
00
(, ), (1.1)() (1.2)y f x y x x X y x a '=<≤??=? 的数值解的方法,="" 次年又把它推广到二阶方程。="" 欧拉的想法如下:我们选择="" 0h="">, 然后在 00x x x h ≤≤+情况下用解函数的切线
0000() () (, ) l x y x x f x y =+-
代替解函数。这样对于点
10x x h =+
就得到
1000(, ) y y hf x y =+。
在 11(, ) x y 重复如上的程序再次计算新的方向就会得到所谓的递推公式:
11, (, ), m m m m m m x x h y y hf x y ++=+=+
这就是 Euler 方法。通过连接所有这些切线得到的函数被称为 Euler 折线。如 果我们令 0h →, 这些折线就会越来越接近解函数。
Euler 折线法是最早出现的, 虽然它亦是常微分方程初值问题的最简单的数 值解法, 但它的一些特性和研究方法对于更复杂的方法却具有普遍意义。几十 年后,法国数学家 A .L .Cauchy(1789.8— 1857.5) 在历史上首次研究了常微分 方程的局部性态。对于给定的初值问题 (1.1)和 (1.1),在 (, ) f x y 连续可微的假 设下, 他用类似于欧拉折线的方法构造逼近解, 利用微分中值定理估计逼近解 之间差的上界,严格证明了以 0x 为中心的一个小区间上逼近解收敛, 其极限函 数即为所提问题的解。同时 Cauchy 指出,这种方法也适用于常微分方程组,所 以欧拉方法有时又称 Euler-Cauchy 折线法。
2.2 Runge-Kutta方法
德国数学家 C.D.T.Runge(1856— 1927) 是数值方法发展史上具有里程碑作 用的人物。 1895年, 他在 Hanover 发表了关于微分方程数值解法的经典论文 《常 微分方程数值解法》 , 此文成为常微分方程 Runge-Kutta 方法的发端。此后, Runge 结合教学活动积极投身于发展一般的数值分析特别是各种实际应用中的 Runge-Kutta 方法 (严格来说, 此方法在 Kutta 做出工作后才能称作 Runge-Kutta 方法) 。
Runge 伟大创造的思想是什么呢?他的灵感来自于初值问题 (1.1) 和 (1.2) 与积分问题
0000() () () x h
x y x h y x f x dx f y ++=+
?(此时 与 无关) (1.3)
之间的对比, 显然, 等式 (1,3) 右侧数值积分的精确度决定 0() y x h +的精确度, Runge 发现, Euler方法采用的是左矩形公式
000() () x h
x f x dx hf x +≈?,
即用高为 0() f x 宽为 h 的矩形代替数值积分, 而这个公式的精确度并不高。 因此 他说:最好通过插入上述 Euler 步骤的结果来代替未出现的 y 值, 把精度更高 的中点法则和梯形法则拓展到微分方程。
10000011:(, (, )) 22
M y y hf x h y hf x y ==+++, 1000000011:((, ) (, (, ))) 22T y y h f x y f x h y hf x y ==++++。 其中 M 和 T 分别表示用中点法和梯形法算得的数值积分。与他的后继者一样, Runge 用 Taylor 展开式表明上述两方法的局部误差是,方法的阶为 2。不过他 的 梦 想 却 是 使 用 具 有 更 高 精 度 的 Simpson 法 则 。 但 是 众 所 周 知 , () /3R M T M =+-的微小变化往往易产生假象,令人误以为可以获得更高的阶。 Taylor 级数展开表明,如果 f 依赖于 y ,事实上这个表达形式只是 2 阶的。接 着 Runge 发现,通过重复使用 Euler 步骤对梯形法则做些许调整,会使形式 () /3R M T M =+-'成为 3 阶方法。 Runge 还把他的方法及方法的展开式拓展到 微分方程组。
1900 年, Runge 的同胞 K .Heun 评论说, Runge 获得的上述方法是归纳性的 而且是令人费解的, 他提出采用 “更具一般性” 的 Gauss 方法。 于是一般的 Gauss 求积公式
1001() s
i i i y y h b f x c h ==++∑,
被扩充为
100202021303032(, ),
(, ), (, ),
k f x y k f x c h y c hk k f x c h y c hk ==++=+
+ 101s
i i i y y h b k ==+∑ (1.4)
把(1.4)式的右端进行二元 Taylor 展开后与 () y x h +的 Taylor 展开式的对应 项的系数比较,适当选取参数使方法具有尽可能高的精度。
Runge 的另一个同胞 W.M.Kutta(1867— 1944) , 1894到 1909年在 Munich 做 助教,在那里受到 Runge 文章的吸引并在 Heun 论文的激励下发表了他的研究结 果。 他认为, 为什么不让已经求得的导数值全部进入到新的求值点的计算中呢?
基于这样的想法, Heun 格式就被 Kutta 代替为如下格式:
100202021130303113224040411422433(, ),
(, ),
(, ),
(, )
k f x y k f x c h y a hk k f x c h y a hk a hk k f x c h y a hk a hk a hk ==++=+++=+++
+
101 s i i i y y h b k ==+∑
(其中 s 称为级 ) 这个格式在满足所需的阶条件上能够允许更多的自由度。 在古典 的 Runge-Kutta 方法中, 对系数的选择极大地取决于由这些系数构成的方法是否 方便进行桌上计算, 而所谓的古典方法是指在前计算机时代得到的方法。 而这些 方法对于在自动计算机上使用则未必是最适合的方法。
显然 Runge-Kutta 方法是一种特殊的单步方法, 事实上这个方法可以看作在 1(, ) m m x x +上取若干条积分曲线的若干个点的切线斜率, 再进行一次 (或多次) 算 术(或加权)平均后产生的新斜率,再按这个斜率从 (, ) m m x y 出发,以直线代曲 线向前推进一步的过程。与 Taylor 展开方法相比, Runge-Kutta 方法不用增加 微商 (, ) f x y 的次数就可以得到较高的阶。
在 Runge 于 1895年提出 Runge-Kutta 方法的雏形的时候,他给出的方法是 2级 2阶和 4级 3阶的; Heun 在 1900年获得的两个方法分别是 3级 3阶和 8级 4阶的;而 Kutta 在 1901年得到的两个方法则分别是 4级 4阶和 6级 5阶的。 因为 Runge 是提出这个方法雏形的人而 Kutta 则是得到了此方法在前计算时代最 高的阶,所以这个方法被称为 Runge-Kutta 方法。到 1925 年,另外一位数学家 Nystr ?m 也得到了 6级 5阶的方法。 此后直到计算机诞生也未产生新的重要成果, 所以 Nystr ?m 在 1925年得到 6级 5阶方法的论文有时也被称为古典 Runge-Kutta 方法与现代 Runge-Kutta 方法的分水岭。
2.3现代 Runge-Kutta 方法
Runge-Kutta 方法在创立之初并未达到完善, 后来又经过大量数学家的共同 努力才日趋成熟。 20世纪 60年代,新西兰数学家 J.C.Butcher (1933—)给出 了现代 Runge-Kutta 方法的一般形式:
1111;
(, ), , 1, 2,
, ; , 1, 2, , . s n n i i i s i n i n ij j j s i ij j y y h b k k f x c h y h a k i j s c a i s +===?=+??
??=++=??
?==???∑∑∑ (1.5)
为了展示 (1.5) 式中出现的系数, Butcher 给出了后来以其名字命名的 Butcher 点阵:
1
111212
212221
212s s s s s ss s
c a a a c a a a c a a a b b
b 我们分别把 s 维向量 c 、 b 以及 s s ?矩阵 A 定义为
[][]1212, , , , , , , , . T s s ij c c c c b b b b A a ??===??
显然, 一 s 级 Runge-Kutta 方法被其 Butcher 点阵完全确定,这样,研究
某一 Runge-Kutta 方法的性质就转化为研究与其相对应的矩阵 A 的性质。 此点阵 为研究 Runge-Kutta 方法的性质提供了方便的途径。在计算机诞生之前, Runge-Kutta 方法被禁锢在只能用手进行计算的实际问题上, 但是随着计算机的 出现, Runge-Kutta 方法呈现出新的史无前例的重要性。能够解决的问题变得越 来越大、 越来越复杂, 自动化的误差检测和步长控制不但变得可行而且变得必要 了, Runge-Kutta 方法的发展不但表现在理论上而且表现在技术上。
除了 Runge-Kutta 方法在微分方程求解中扮演的传统角色外, 人们发现相关
类型的初值问题可以用 Runge-Kutta 方法或适合于更一般问题的 Runge-Kutta 方法求解,比如 Runge-Kutta 方法被应用到了 Hamilton 系统中。
二、实际应用 [2]
2.1 经典四阶 Runge-Kutta 方法
经典的四阶 Runge-Kutta 方法是:
112341213243(22) 6(, )
11(, ) 2211(, ) 22(, )
n n n n n n n n n n h y y k k k k k f x y k f x h y hk k f x h y hk k f x h y hk +?=++++??=???=++???=++??=++?? 它 的 局 部 截 断 误 差 51() n T O h +=, 所 以 为 四 阶 方 法 , 这 是 最 常 用 的 四 阶 Runge-Kutta 方法,数学库中都有用此方法求解初值问题的软件,这种方法的优 点是精度较高,缺点是每步要计算四个函数值,计算量较大。
2.2 两种群共存模型的 Runge-Kutta 法解答
2.2.1 问题提出
在自然界中处于同一环境下有多个种群生存,它们之间存在着一定的关系, 它们之间可能是食饵与捕食者, 或是相互依存, 或是相互竞争或其它的关系。 简 单的食饵与捕食者模型是由两位数学生态学的先驱者 A.J.Lotka 与意大利数学 家 V.Volterra 各自建立, 因此模型被命名为 Lotka-Volterra 模型。 在该文献中, 在 Lotka-Volterra 模型的基础上建立考虑捕食者另有食物来源,并且把被捕食 种群作为生存资源, 限制捕食种群增长的两种群共存的数学模型, 并对模型进行 分析和模型解释。两种群相互共生存的现象是普遍存在的 , 研究其种群的数量演 变规律具有重要意义, 它有助于我们对种群变化情况的了解和正确解释, 可以帮 助我们合理解决这些生态问题。 下面是该文献的两种群相互共存的微分方程模型。
2.2.2 模型建立与分析
①假设有种群 A 和种群 B ,种群 A 是种群 B 的食物来源,但不是唯一的食物 来源,种群 B 还另有食物来源。
②假设 1M 是环境容许的种群 A 的最大数量, 2M 是环境容许的种群 B 的最大 数量。
③假设种群 A 能独立生存,设其固有增长率为 1μ,由于种群 A 是种群 B 的 食物来源,故设 10λ>为单位数量的种群 B (相对 2M 而言)掠取 1λ倍的单位种群 B 的量(相对 1M 而言) ,还考虑到种群 A 的自身的阻滞作用,因此设:
21111121() 1y y y t y M M μλ??=-- ??
?, 1() y t 表示在 t 时刻种群 A 的数量。 ④因为种群 B 另有食物来源,设其增长率为 2μ,假设种群 A 作为生存资源 限制种群 B 的增长,设对其增长率的影响率为 221
y y λ-, (20λ>) (表示种群 A 对 种群 B 的影响,若种群 A 充足, 21
y y 小,则对种群 B 的增长影响较小;若种群 A 不充足, 21
y y 大,则对种群 B 的增长影响较大) ,不考虑种群 B 自身内部竞争,因 此设 222221() 1y y t y y μλ??=- ??
?, 2() y t 表示在 t 时刻种群 B 的数量。 根据上面的假设可以得到相互共存的两种群的微分方程模型。
1211112122222111dy y y y dt M M dy y y dt y μλμλ???=--? ????????=- ?????
(2.1) 其初始条件为 110(0)y y =, 220(0)y y =。
将(2.1)可归结为如下形式的一阶常微分方程组
11122212(, , ) (, , ) dy F t y y dt dy F t y y dt
?=????=?? (2.2) 因为此方程组没有解析解, 故采用标准的四阶 Runge-Kutta 方法求出数值解, 求 解公式如下:
1, 11, 112131412, 12, 1222324211, 2, 21, 112, 1231, 212, 2241, 312, (22) 6(22) 6(, , )
111, , 222111, , 22211, , 22n n n n j j n n n j j n n n j j n n n j j n n n h y y k k k k h y y k k k k k F t y y k F t h y hk y hk k F t h y hk y hk k F t h y hk y ++=++++=++++=??=+++ ?????
=+++ ???=+++()321, 2j hk ????????=????????? ?????
因为 12111111
111212222222222222112121, 2, F y F y y M M y M F y F y y y y y μλμλμλμλμλμ????=--=- ????
???==-??
存在且连续,由定理可知采用标准的四阶 Runge-Kutta 方法是收敛的。
为了对微分方程模型进行数值计算,假设初始值为
12(0)100, (0)200y y ==,
方程组 (2.1) 中 10. 6μ=, 20.4μ=, 10.5λ=, 20.2λ=, 11000M =, 23000M =, 用 Matlab 编程如下:
function dq=myfile_2(t,y)
dq=[0.5*y(1)*(l-0.5*y(2)/3000-y(l)/1000);0.4*y(2)*(l-0.2*y(2)/y(l))];
ode45(@myfile_2,[0,100],[100;200]),grid
于是得 1() y t 、 2() y t 的图形 1如下:
为了体现模型中所设参数 1μ, 2μ, 1λ, 2λ在生态学上的意义,我们再设一 组参数 10.5μ=, 20.3μ=, 10.2λ=, 21.5λ=,用 Matlab 编程如下: function dq=myfile_3(t,y)
dq=[0.5*y(1)*(l-0.2*y(2)/3000-y(l)/1000);0.3*y(2)*(l-1.5*y(2)/y(l))];
ode45(@myfile_3,[0,100],[100:200]),grid
于是得图形 2:
比较上述两个图形可以看出种群 A 与种群 B 的相互影响, 种群 B 的数量与对 其增长率的影响为 221
y y λ-, ()20λ>有关,系数 2λ可以反映种群 A 对种群 B 的影 响程度,系数 2λ越小,种群 A 对种群 B 的增长影响较小,系数 2λ越大,种群 A 对种群 B 的增长影响较大, 种群 B 的数量随着种群 A 的数量的变化而出现不同的 情况变化。系数 1λ可以反映种群 B 对种群 A 的影响程度,系数 1λ越小,种群 B 对种群 A 的增长影响较小, 系数 1λ越大, 种群 B 对种群 A 的增长影响较大, 从而 种群 A 的数量变化出现不同的情况。 从上述图形分析比较可以知道, 通过改变种 群 A 的数量可以使种群 B 的数量发生变化,但最终两种群趋向于稳
定而共存。
三、总结
通过在网上查阅资料, 本文介绍了常微分方程初值解法—— Runge-Kutta 法 的发展历程和其中的主要人物,并将其方法应用到求解建立的常微分方程模型, 通过数值计算作出图形, 并对模型作出分析, 从而帮助我们分析解决实际问题和 实际现象的解释, 有利于我们采取实际有效的措施解决问题。
由于我们在解决现
实中的许多问题往往通过建立数学模型来解决, 而建立的模型很多是微分方程模 型, 并且很难或无法求出解析解, 因此, 运用数值解法对于求解常微分方程模型 的数值解有着重要的意义, 研究微分方程的数值求解具有很大的潜力, 是值得我 们去探讨和研究的。
参考文献:
[1]林立军, 常微分方程数值解法—— Runge_Kutta法的历史浅析 [J], 辽宁师范大学学报 (自 然科学版) , 2003
[2]秦军, Rung-Kutta 法在求解微分方程模型中的应用 [J], 2010
[3]李庆扬 , 王能超 , 易大义 . 数值分析 [M],清华大学出版社, 2001
范文三：数值分析计算方法
《计算方法》实验内容
10000
一．实验一：用两种不同的顺序计算∑n n =1-2≈1. 644834，分析其误差的变化。
1. 实验目的：通过正序反序两种不同的顺序求和，比较不同算法的误差；了解在计算机中大数吃小数的现象，以后尽量避免；体会单精度和双精度数据的差别。
2. 算法描述：累加和s=0；
正序求和：
对于n=1,2,3,......,10000
s+=1.0/(n*n);
反序求和：
对于n=10000,9999,9998,.....,1
s+=1.0/(n*n);
3. 源程序：
#双精度型#
#include void main() { double s=0; int n; for(n=1;n<> s+=1.0/(n*n); printf("正序求和结果是：%lf\n",s); s=0; for(n=10000;n>=1;n--) s+=1.0/(n*n); printf("反序求和结果是:%lf\n",s); } #单精度型# #include void main() { float s=0; int n; for(n=1;n<> s+=1.0/(n*n); printf("正序求和结果是：%f\n",s); s=0; for(n=10000;n>=1;n--) s+=1.0/(n*n); printf("反序求和结果是:%f\n",s); } 4. 运行结果： 双精度型运行结果： 单精度型运行结果： 5. 对算法的理解与分析：舍入误差在计算机中会引起熟知的不稳定，算法不同，肯结果也会不同，因此选取稳定的算法很重要。选取双精度型数据正反序求和时结果一致，但选用单精度型数据时，求和结果不一致，明显正序求和结果有误差，所以第一个算法较为稳定可靠。 二．实验二： 1、拉格朗日插值 作二次插值，并求x 1=-2，x 2=0，x 3=2.75时的函数近似值 2牛顿插值 作五次插值，并求x 1=0.46，x 2=0.55，x 3=0.60时的函数近似值. 1. 实验目的：通过拉格朗日插值和牛顿插值的实例，了解两种求解方法，并分析各自的优缺点。 2. 算法描述： 3. 源程序： 拉格朗日插值： #include #define k 2 void main() { double L,Y,a; double x[3]; double y[3]; for(int p=0;p<> scanf("%lf%lf",&x[p],&y[p]); for(int q=0;q<> { Y=0; scanf("%lf",&a); for(int j=0;j<> { L=1; for(int i=0;i<> if(i!=j) L*=(a-x[i])/(x[j]-x[i]); Y+=L*y[j]; } printf("x=%lf的函数值是:%lf\n",a,Y); } } 牛顿插值： #include void main() { int i,j; double a[6][6],y,s,x; for(j=0;j<> for(i=0;i<> scanf("%lf",&a[i][j]); for(j=2;j<> for(i=j-1;i<> a[i][j]=(a[i][j-1]-a[i-1][j-1])/(a[i][0]-a[i-j+1][0]); for(int k=0;k<> { scanf("%lf",&x); y=0; for(i=0;i<> { s=a[i][i+1]; for(j=i-1;j>=0;j--) s*=x-a[j][0]; y+=s; } printf("结果是：%lf\n",y); } 4. 运行结果： 拉格朗日运行结果： 牛顿插值运行结果： 5. 对算法的理解与分析： （1）拉格朗日插值： 该公式是对一系列点加权求和。内插通常优于外推。选择的区间包括x ，插值效果越好，高次插值通常优于低次插值，但并不是意味着插值次数越高越好。 优点：编程容易实现。 缺点：如果发现当前的插值方法 不够精确，就要增加插值节点的个数。拉格朗日基函数每一个都要重新计算，效率低。 （2）牛顿插值： 优点：如果当前插值方法不稳定，需要增加插值节点个数，只需要计算新家节点所增加的差商即可，之前的计算结果可以被复用，比拉格朗日插值节省计算量，效率高，精度高，更为稳定。 三．实验三：分别用复化梯形公式和复化辛卜生公式计算f(x)=sin(x)/x的积分， 并与准确值比较判断精度。 1. 实验目的：通过实例体会各种算法的精度。熟练掌握复化梯形，复化辛普森，复化柯特斯求积方法的程序。 2. 算法描述： 复化梯形： fc(x)=sinx/x,N等分,a,b 是积分区间; 令s=fc(a)+fc(b), 对于i=a+1.0/N,a+2.0/N,.....(i<> s+=2*fc(i); s*=double(b-a)/(2*N); 输出s; 复化辛普森： fc(x)=sinx/x,N等分,a,b 是积分区间; 令s=fc(a)+fc(b), 对于i=a+1.0/N,a+2.0/N,.....(i<> s+=2*fc(i); 对于i=a+0.5/N,a+1.5/N,.....(i<> s+=4*fc(i); s*=double(b-a)/(6*N); 输出s ； 复化柯特斯： fc(x)=sinx/x,N等分,a,b 是积分区间; 令s=fc(a)+fc(b), 对于i=a+1.0/(4*N),a+2.0/(4*N),.....(i<> s+=32*fc(i); 对于i=a+2.0/(4*N),a+6.0/(4*N),.....(i<> s+=12*fc(i); 对于i=a+1.0/N,a+2.0/N,.....(i<> s+=14*fc(i); s*=double(b-a)/(90*N); 输出s; 3. 源程序： #include #include #define N 8 #define a 0 #define b 1 double fc(double x) { double y; y=sin(x)/x; return y; } void tx() double s,i; s=fc(a)+fc(b); for(i=a+1.0/N;i s+=2*fc(i); s*=double(b-a)/(2*N); printf("复合梯形8等分求积结果是：%lf\n",s); } void xps() { double s,i; s=fc(a)+fc(b); for(i=a+1.0/N;i s+=2*fc(i); for(i=a+0.5/N;i s+=4*fc(i); s*=double(b-a)/(6*N); printf("复合辛普森8等分求积结果是：%lf\n",s); } void kts() { double s,i; s=fc(a)+fc(b); for(i=a+1.0/(4*N);i s+=32*fc(i); for(i=a+2.0/(4*N);i s+=12*fc(i); for(i=a+1.0/N;i s+=14*fc(i); s*=double(b-a)/(90*N); printf("复合柯特思斯8等分求积结果是：%lf\n",s); } void main() { char d; scanf("%c",&d); if(d=='t') tx(); if(d=='x') xps(); if(d=='k') kts(); } 4. 运行结果： 5. 对算法的理解与分析：他们都适用于求数值积分，而且都能提高计算的精度，他们都是在每个小区间上在使用梯形，辛普森，柯特斯求积公式，其中复化辛普森和复化柯特斯的收敛速度较快。 四．实验四：用改进欧拉方法解初值问题y’=x+y; y(0)=1。0<><1,取步长h=0.1计算，并与准确值 y="-x-1-2ex"> 1. 实验目的：熟悉常微分方程初值问题的求解方法；熟悉其算法； 2. 算法描述： 3. 源程序： 4. 运行结果： 5. 对算法的理解与分析： 五．实验五：分别用下列方法求f(x)=x3-3x-1=0在x 0=2附近的根。根的准确值为x*=1.87938524…，要求准确到四位有效数字，并对比各种算法的计算量。 （1）二分法；（2）简单迭代法；（3）牛顿迭代法 1. 实验目的：通过对二分法，简单迭代法和牛顿迭代法的编程和上机实验，进一步体会他们在方程求根中的不同特点；比较三者的计算速度和计算精度。 2. 算法描述： 二分法： 给定区间[a,b],e=pow(10,-4)*0.5 (1)令c=(a+b)/2,保持b-a<> (2)如果f(c)*f(a)=0,输出c; 如果f(c)*f(a)<> 否则a=c,执行（1） 简单迭代法： 3. 源程序： 二分法： #include #include double fc(double x) { double y; y=x*x*x-3*x-1; return y; } void main(void) { double fa=fc(1),fb=fc(3),a=1,b=3,f,x0; int k=0; for(;b>a&&b-a>=pow(10,-4)*0.5;) { f=fc((a+b)/2); if(f==0) { x0=(a+b)/2; break; } else if(fa*f<> b=(a+b)/2; else a=(a+b)/2; k++; } x0=(a+b)/2; printf("近似根是：%lf 准确根是：1.87938524\n",x0); printf("迭代次数是：%d\n",k); } 简单迭代法： #include #include double Iterate1(double x) { double y; y=pow((3*x+1),1.0/3); return y; } double Derivative1(double x) { double y; y=pow((3*x+1),-2.0/3); return y; } double Iterate2(double x) { double y; y=(1-x*x*x)/3.0; return y; } double Derivative2(double x) { double y; y=-x*x; return y; } double Iterate3(double x) { double y; y=(3*x+1)/(x*x); return y; } double Derivative3(double x) { double y; y=(-3*x-2)/(x*x*x); return y; } void main(void) { double x2=2.0,x1; int k=0; double a1=Derivative1(x2); if(fabs(a1)<> { printf("迭代公式是x=pow((3*x+1),1.0/3);\n"); do { x1=x2; x2=Iterate1(x1); k++; }while(fabs(x2-x1)>=pow(10,-4)*0.5); printf("近似根是：%lf 准确根是：1.87938524\n",x1); printf("迭代次数是：%d\n",k); } double a2=Derivative2(x2); if(fabs(a2)<> { printf("迭代公式是x=(1-x*x*x)/3.0;\n"); do { x1=x2; x2=Iterate2(x1); k++; }while(fabs(x2-x1)>=pow(10,-4)*0.5); printf("近似根是：%lf 准确根是：1.87938524\n",x1); printf("迭代次数是：%d\n",k); } double a3=Derivative3(x2); if(fabs(a3)<> { printf("迭代公式是x=(3*x+1)/(x*x);\n"); do { x1=x2; x2=Iterate3(x1); k++; }while(fabs(x2-x1)>=pow(10,-4)*0.5); printf("近似根是：%lf 准确根是：1.87938524\n",x1); printf("迭代次数是：%d\n",k); } } 牛顿迭代法： #include #include double Iterate(double x) { double y; y=x*x*x-3*x-1; return y; } double Derivative(double x) { double y; y=3*x*x-3; return y; } void main() { double x0=2.0,x1; double f0,f01,f1,f11; f0=Iterate(x0); f01=Derivative(x0); int k=0; x1=x0-f0/f01; for(;fabs(x0-x1)>=pow(10,-4)*0.5;f0=f1,f01=f11) { x0=x1; x1=x0-f0/f01; f1=Iterate(x1); f11=Derivative(x1); k++; if(f11==0) { printf("牛顿迭代法失效\n"); break; } } printf("近似根是：%lf 准确根是：1.87938524\n",x1); printf("迭代次数是：%d\n",k); } 4. 运行结果： 二分法运行结果： 简单迭代法运行结果： 牛顿迭代法运行结果： 插值法 1. 下列数据点的插值 x 0 1 4 9 16 25 36 49 64 y 0 1 2 3 4 5 6 7 8 可以得到平方根函数的近似，在区间[0,64]上作图. (1)用这9个点作8次多项式插值Ls(x). (2)用三次样条(第一边界条件)程序求S(x). 从得到结果看在[0,64]上，哪个插值更精确;在区间[0,1]上，两种插值哪个更精确, 解:(1)拉格朗日插值多项式，求解程序如下 syms x l; x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; n=length(x1); Ls=sym(0); for i=1:n l=sym(y1(i)); for k=1:i-1 l=l*(x-x1(k))/(x1(i)-x1(k)); end for k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k)); end Ls=Ls+l; end Ls=simplify(Ls) %为所求插值多项式Ls(x). 输出结果为 Ls = -24221063/63504000*x^2+95549/72072*x-1/3048192000*x^8-2168879/435456000 *x^4+19/283046400*x^7+657859/10886400*x^3+33983/152409600*x^5-13003/2395008 000*x^6 (2)三次样条插值，程序如下 x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; x2=[0:1:64]; y3=spline(x1,y1,x2); p=polyfit(x2,y3,3); %得到三次样条拟合函数 S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x) 输出结果为: S = 2288075067923491/73786976294838206464-2399112304472833/5764607523034234 88*x+4552380473376713/18014398509481984*x^2+999337332656867/112589990684262 4*x^3 (3)在区间[0,64]上，分别对这两种插值和标准函数作图， plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y') 蓝色曲线为y=?x函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线 100 80 60 40 20 0 -20010203040506070 可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。 在[0,1]区间上由上图看不出差别，不妨代入几组数据进行比 较 ，取x4=[0:0.2:1] x4=[0:0.2:1]; sqrt(x4) %准确值 subs(Ls,'x',x4) %拉格朗日插值 spline(x1,y1,x4) %三次样条插值 运行结果为 ans = 0 0.4472 0.6325 0.7746 0.8944 1.0000 ans = 0 0.2504 0.4730 0.6706 0.8455 1.0000 ans = 0 0.2429 0.4630 0.6617 0.8403 1.0000 从这几组数值上可以看出在[0,1]区间上，拉格朗日插值更精确。 数据拟合和最佳平方逼近 2. 有实验给出数据表 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次、4次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。 解:(1)三次拟合，程序如下 sym x; x0=[0.0 0.1 0.2 0.3 0.5 0.8 1.0]; y0=[1.0 0.41 0.50 0.61 0.91 2.02 2.46]; cc=polyfit(x0,y0,3); S3=cc(1)+cc(2)*x+cc(3)*x^2+cc(4)*x^3 %三次拟合多项式 xx=x0(1):0.1:x0(length(x0)); yy=polyval(cc,xx); plot(xx,yy,'--'); hold on; plot(x0,y0,'x'); xlabel('x'); ylabel('y'); 运行结果 S3 = -7455778416425083/1125899906842624+1803512222945437/140737488355328*x-655705 280524945/140737488355328*x^2+4172976892199509/4503599627370496*x^3 图像如下 2.5 2 1.5 y 1 0.5 000.10.20.30.40.50.60.70.80.91 x (2)4次多项式拟合 sym x; x0=[0.0 0.1 0.2 0.3 0.5 0.8 1.0]; y0=[1.0 0.41 0.50 0.61 0.91 2.02 2.46]; cc=polyfit(x0,y0,4); S3=cc(1)+cc(2)*x+cc(3)*x^2+cc(4)*x^3+cc(5)*x^4 xx=x0(1):0.1:x0(length(x0)); yy=polyval(cc,xx); plot(xx,yy,'r'); hold on; plot(x0,y0,'x'); xlabel('x'); ylabel('y'); 运行结果 S3 = 3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931 990070637/281474976710656*x^2-5965836931688425/1125899906842624*x^3+84912754646 50307/9007199254740992*x^4 图像如下 2.5 2 1.5 y 1 0.5 000.10.20.30.40.50.60.70.80.91 x (3)另一个拟合曲线, 新建一个M-file 程序如下: 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 在命令窗口中输入以下的命令: 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 on; plot(xx,yy,'b',x,y,'.'); 图像如下 2.5 2 1.5 y 1 0.5 000.10.20.30.40.50.60.70.80.91 x 解线性方程组的直接解法 3. 线性方程组Ax=b的A及b为 10 7 8 7321 7 5 6 5231A=,,，b=,,，则解x=,,.用MATLAB内部函数求detA8 6 10 9331 7 5 9 10311 及A的所有特征值和cond(A)2.若令 10 7 8.1 7.2 7.08 5.04 6 5A+δA=,,，求解(A+δA)(x+δx)=b,输出向量x和8 5.98 9.89 9 6.99 5 9 9.98 ||δx||2.从理论结果和实际计算两方面分析线性方程组Ax=b解 得相对误差||δx||2/||x||2及A的相对误差||δA||2/||A||2的 关系. 解:(1)程序如下 clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; det(A) cond(A,2) eig(A) 输出结果为 ans = 1 ans = 2.9841e+003 ans = 0.0102 0.8431 3.8581 30.2887 (2)程序如下 A=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98]; b=[32 23 33 31]'; x0=[1 1 1 1]'; x=A\b %扰动后方程组的解 x1=x-x0 %δx的值 norm(x1,2) %δx的2-范数 运行结果为 x = -9.5863 18.3741 -3.2258 3.5240 x1 = -10.5863 17.3741 -4.2258 2.5240 ans = 20.9322 程序如下 A=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98]; A0=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; b=[32 23 33 31]'; x0=[1 1 1 1]'; x=A\b; x1=x-x0; norm(x1,2); A1=A-A0 ; %δA的值 norm(x1,2)/norm(x0,2) % ||δx||2/||x||2的值 |A||2的值 norm(A1,2)/norm(A0,2) %||δA||2/|输出结果为 ans = 10.4661 ans = 0.0076 可见A相对误差只为0.0076，而得到的结果x的相对误差就达到了 10.4661，该方程组是病态的，A的条件数为2984.1远远大于1，当 A只有很小的误差就会给结果带来很大的影响。 非线性方程数值解法 4. 求下列方程的实根 (1) x^2-3x+2-e^x=0; (2) x^3+2x^2+10x-20=0. 要求:(1)设计一种不动点迭代法，要使迭代序列收敛，然后再用 斯特芬森加速迭代，计算到|x(k)-x(k-1)|<> 牛顿迭代，同样计算到|x(k)-x(k-1)|<> 各次迭代值和迭代次数k，比较方法的优劣。 解:(1)先用画图的方法估计根的范围 ezplot('x^2-3*x+2-exp(x)'); grid on; 2-3 x+2-exp(x)x 50 0 -50 -100 -150 -200 -6-4-20246 x 可以估计到方程的根在区间(0,1);选取迭代初值为x0=0.5; 构造不动点迭代公式x(k+1)=( x(k)^2+2-e^x(k))/3; ψ(x)= ( x^2+2-e^x)/3; 当0<><><><1;><> 程序如下: format long; f=inline('(x^2+2-exp(x))/3') disp('x='); x=feval(f,0.5); disp(x); Eps=1E-8; i=1; while 1 x0=x; i=i+1; x=feval(f,x); disp(x); if abs(x-x0) break; end end i,x 运行结果为 f = Inline function: f(x) = (x^2+2-exp(x))/3 x= 0.20042624309996 0.27274906509837 0.25360715658413 0.25855037626494 0.25726563633509 0.25759898516219 0.25751245451483 0.25753491361525 0.25752908416796 0.25753059723833 0.25753020451046 0.25753030644564 0.25753027998767 0.25753028685501 i = 14 x = 0.25753028685501 斯特芬森加速法，程序如下: format long; f=inline('x-((x^2+2-exp(x))/3-x)^2/((((x^2+2-exp(x))/3)^2+2-exp((x^2+2- exp(x))/3))/3-2*(x^2+2-exp(x))/3+x)'); disp('x='); x=feval(f,0.5); disp(x); Eps=1E-8; i=1; while 1 x0=x; i=i+1; x=feval(f,x); disp(x); if abs(x-x0) break; end end i,x 运行结果为x= 0.25868442756579 0.25753031771981 0.25753028543986 0.25753028543986 i = 4 x = 0.25753028543986 牛顿迭代法，程序如下: format long; x=sym('x'); f=sym('x^2-3*x+2-exp(x)'); df=diff(f,x); FX=x-f/df; Fx=inline(FX); disp('x='); x1=0.5; disp(x1); -8; Eps=1E i=0; while 1 x0=x1; i=i+1; x1=feval(Fx,x1); disp(x1); if abs(x1-x0) break; end end i,x1 运行结果如下: x= 0.50000000000000 0.25368870241829 0.25752890079471 0.25753028543968 0.25753028543986 i = 4 x1 = 0.25753028543986 (2) 先用画图的方法估计根的范围 ezplot('x^3+2*x^2+10*x-20'); grid on; 32+2 x+10 x-20x 400 300 200 100 0 -100 -200 -6-4-20246 x 根大约在区间(1,2);选取初值x0=1.5; 构造不动点迭代公式x(k+1)=(-2x(k)^2-10x(k)+20)^1/3; ψ(x)=(-2x^2-10x+20)^1/3; 程序如下: format long; f=inline('(-2*x^2-10*x+20)^1/3') disp('x='); x=feval(f,1.5); disp(x) Eps=1E-8; i=1; while 1 x0=x; i=i+1; x=feval(f,x); disp(x); if abs(x-x0)>1E10 break; end if abs(x-x0) break; end end i,x 运行结果: f = Inline function: f(x) = (-2*x^2-10*x+20)^1/3 x= 0.16666666666667 6.09259259259259 -38.38843164151806 -8.478196837919431e+002 -4.763660785374071e+005 -1.512815059604763e+011 i = 6 x = -1.512815059604763e+011 迭代6次后x的值大得令人吃惊,表明构造的式子并不收敛. 也无法构造出收敛的不动点公式 牛顿迭代法，程序如下: format long; x=sym('x'); f=sym('x^3+2*x^2+10*x-20'); df=diff(f,x); FX=x-f/df; Fx=inline(FX); disp('x='); x1=0.5; disp(x1); Eps=1E-8; i=0; while 1 x0=x1; i=i+1; x1=feval(Fx,x1); disp(x1); if abs(x1-x0) break; end end i,x1 运行结果: x= 1.50000000000000 1.37362637362637 1.36881481962396 1.36880810783441 1.36880810782137 i = 4 x1 = 1.36880810782137 比较三种方法，牛顿法的收敛性比较好，相比不动点迭代法要构造出 收敛的公式比较难，牛顿法迭代次数也较少，收敛速度快，只是对初 值的要求很高，几种方法各有利弊，具体采用哪种也需因题而异。 常微分方程初值问题数值解法 5. 给定初值问题 y’=-50y+50x^2+2x,0?x?1; y(0)=1/3; 用经典的四阶R-K方法解该问题，步长分别取h=0.1,0.025,0.01, 计算并打印x=0.1i(i=0,1,…,10)各点的值，与准确值 y(x)=1/3e^(-50x)+x^2比较。 解:取步长h=0.1，程序如下: %经典的四阶R-K方法 clear; F='-50*y+50*x^2+2*x'; a=0;b=1; h=0.1; n=(b-a)/0.1; X=a:0.1:b; Y=zeros(1,n+1); Y(1)=1/3; for i=1:n x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6; end %准确值 temp=[]; f=dsolve('Dy=-50*y+50*x^2+2*x','y(0)=1/3','x'); df=zeros(1,n+1); for i=1:n+1 temp=subs(f,'x',X(i)); df(i)=double(vpa(temp)); end disp(' 步长 四阶经典R-K法 准确值'); disp([X',Y',df']); 运行结果: 步长 四阶经典R-K法 准确值 1.0e+010 * 0 0.00000000003333 0.00000000003333 0.00000000001000 0.00000000046055 0.00000000000122 0.00000000002000 0.00000000630625 0.00000000000400 0.00000000003000 0.00000008640494 0.00000000000900 0.00000000004000 0.00000118436300 0.00000000001600 0.00000000005000 0.00001623545110 0.00000000002500 0.00000000006000 0.00022256067134 0.00000000003600 0.00000000007000 0.00305093542778 0.00000000004900 0.00000000008000 0.04182323921740 0.00000000006400 0.00000000009000 0.57332690347809 0.00000000008100 0.00000000010000 7.85935630083771 0.00000000010000 %画图观察结果 figure; plot(X,df,'k*',X,Y,'--r'); grid on; title('四阶经典R-K法解常微分方程'); legend('准确值','四阶经典R-K法'); 10四阶经典R-K法解常微分方程x 108 准确值 四阶经典R-K法7 6 5 4 3 2 1 000.10.20.30.40.50.60.70.80.91 当x值接近1的时候，偏离准确值太大。 当步长h=0.025时，将上面程序中的h改为0.025即可，运行结果: 步长 四阶经典R-K法 准确值 0 0.33333333333333 0.33333333333333 0.10000000000000 0.10313524034288 0.01224598233303 0.20000000000000 0.04428527327599 0.04001513330992 0.30000000000000 0.05196795755507 0.09000010196744 0.40000000000000 0.09395731149439 0.16000000068705 0.50000000000000 0.16034531435757 0.25000000000463 0.60000000000000 0.24808570130557 0.36000000000003 0.70000000000000 0.35624188472758 0.49000000000000 0.80000000000000 0.48452590661627 0.64000000000000 0.90000000000000 0.63284923300751 0.81000000000000 1.00000000000000 0.80118464374206 1.00000000000000 图像如下: 四阶经典R-K法解常微分方程 1 准确值 0.9四阶经典R-K法 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 000.10.20.30.40.50.60.70.80.91 当步长h=0.01时，将上面程序中的h改为0.01即可，运行结果: 步长 四阶经典R-K法 准确值 0 0.33333333333333 0.33333333333333 0.10000000000000 0.20235720486111 0.01224598233303 0.20000000000000 0.12881700190791 0.04001513330992 0.30000000000000 0.09799182667850 0.09000010196744 0.40000000000000 0.10094946775024 0.16000000068705 0.50000000000000 0.13227011975470 0.25000000000463 0.60000000000000 0.18866520287199 0.36000000000003 0.70000000000000 0.26813930278431 0.49000000000000 0.80000000000000 0.36948166028319 0.64000000000000 0.90000000000000 0.49195762199475 0.81000000000000 1.00000000000000 0.63512142167910 1.00000000000000 四阶经典R-K法解常微分方程 1 准确值 0.9四阶经典R-K法 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 000.10.20.30.40.50.60.70.80.91 《数值分析》计算实习题 姓名: 学号: 班级: 第二章 1、程序代码 Clear;clc; x1=[0.2 0.4 0.6 0.8 1.0]; y1=[0.98 0.92 0.81 0.64 0.38]; n=length(y1); c=y1(:); for j=2:n %求差商 for i=n:-1:j c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1)); end end syms x df d; df(1)=1;d(1)=y1(1); for i=2:n %求牛顿差值多项式 -x1(i-1)); df(i)=df(i-1)*(x d(i)=c(i-1)*df(i); end P4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数 pp=csape(x1,y1, 'variational');%调用三次样条函数 q=pp.coefs; q1=q(1,:)*[(x-.2)^3;(x-.2)^2;(x-.2);1]; q1=vpa(collect(q1),5) q2=q(1,:)*[(x-.4)^3;(x-.4)^2;(x-.4);1]; q2=vpa(collect(q2),5) q3=q(1,:)*[(x-.6)^3;(x-.6)^2;(x-.6);1]; q3=vpa(collect(q3),5) q4=q(1,:)*[(x-.8)^3;(x-.8)^2;(x-.8);1]; q4=vpa(collect(q4),5)%求解并化简多项式 2、运行结果 P4 = 0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x - 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784 q1 = - 1.3393*x^3 + 0.80357*x^2 - 0.40714*x + 1.04 q2 = - 1.3393*x^3 + 1.6071*x^2 - 0.88929*x + 1.1643 q3 = - 1.3393*x^3 + 2.4107*x^2 - 1.6929*x + 1.4171 q4 = - 1.3393*x^3 + 3.2143*x^2 - 2.8179*x + 1.8629 1.1 已知的点 1牛顿插值多项式 三次样条函数 0.9 0.8 0.7 0.6 0.5 0.4 00.10.20.30.40.50.60.70.80.91 3、问题结果 Px()4次牛顿差值多项式= 0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x 4 - 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784 三次样条差值多项式 Qx(), 32,,，,，,1.33930.803570.407141.04,[0.2,0.4]xxxx ,32,，,，,1.33931.60710.889291.1643,[0.4,0.6]xxxx,,32,，,，,1.33932.41071.69291.4171,[0.6,0.8]xxxx,32,,，,，,1.33933.21432.81791.8629,[0.8,1.0]xxxx, 第三章 1、程序代码 Clear;clc; x=[0 0.1 0.2 0.3 0.5 0.8 1]; y=[1 0.41 0.5 0.61 0.91 2.02 2.46]; p1=polyfit(x,y,3)%三次多项式拟合 p2=polyfit(x,y,4)%四次多项式拟合 y1=polyval(p1,x); y2=polyval(p2,x);%多项式求值 plot(x,y,'c--',x,y1,'r:',x,y2,'y-.') p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。 y3=polyval(p3,x); plot(x,y,'c--',x,y1,'r:',x,y2,'y-.',x,y3,'k--')%画出四种拟合曲线 2、运行结果 p1 = -6.6221 12.8147 -4.6591 0.9266 p2 = 2.8853 -12.3348 16.2747 -5.2987 0.9427 p3 = 3.1316 -1.2400 0.7356 3 数据曲线 三次多项式拟合 2.5四次多项式拟合 二次多项式拟合 2 1.5 1 0.5 0 00.10.20.30.40.50.60.70.80.91 3、问题结果 32三次多项式拟合P1= -6.622112.81474.65910.9266xxx，,， 432四次多项式拟合P2= 2.885312.334816.27475.29870.9427xxxx,，,， 2二次多项式拟合P3= 3.13161.24000.7356xx,， 第四章 1、程序代码 1)建立函数文件f.m: function y=fun(x); y=sqrt(x)*log(x); 2)编写程序: a. 利用复化梯形公式及复化辛普森公式求解: Clear;clc; h=0.001;%h为步长,可分别令h=1,0.1,0.01,0.001 n=1/h;t=0;s1=0;s2=0; for i=1:n-1 t=t+f(i*h); end T=h/2*(0+2*t+f(1));T=vpa(T,7) %梯形公式 for i=0:n-1 s1=s1+f(h/2+i*h); end for i=1:n-1 s2=s2+f(i*h); end S=h/6*(0+4*s1+2*s2+f(1));S=vpa(S,7) %辛普森公式 a’复化梯形公式和复化辛普生公式程序代码也可以是: Clear;clc; x=0:0.001:1; %h为步长,可分别令h=1,0.1,0.01,0.001 y=sqrt(x).*log(x+eps); T=trapz(x,y); T=vpa(T,7) (只是h=1的运行结果不一样，T=1.110223*10^(-16)，而其余情况结果都相同) Clear;clc; f=inline('sqrt(x).*log(x)',x); S=quadl(f,0,1); S=vpa(S,7) b. 利用龙贝格公式求解: Clear;clc; m=14;%m+1即为二分次数 h=2; for i=1:m h=h/2;n=1/h;t=0; for j=1:n-1 t=t+f(j*h); end T(i)=h/2*(0+2*t+f(1));%梯形公式 end for i=1:m-1 for j=m:i+1 T(j)=4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1); %通过不断的迭代求得T(j)，即T表的对角线上的元素。 end end vpa(T(m),7) 2、运行结果 T = -0.4443875 S = -0.4444345 ans = -0.4444414 3、问题结果 a. 利用复合梯形公式及复合辛普森公式求解: 步长h 1 0.1 0.01 0.001 梯形求积0 -0.4170628 -0.4431179 -0.4443875 T= [1.110223*10^(-16)] 辛普森求-0.3267527 -0.4386308 -0.4441945 -0.4444345 积S= b. 利用龙贝格公式求解: 通过15次二分，得到结果:-0.4444414 第五章 1、程序代码 (1)LU分解解线性方程组: Clear;clc; A=[10 -7 0 1 -3 2.099999 6 2 5 -1 5 -1 2 1 0 2]; b=[8;5.900001;5;1]; [m,n]=size(A); L=eye(n); U=zeros(n); flag='ok'; for i=1:n U(1,i)=A(1,i); end for r=2:n L(r,1)=A(r,1)/U(1,1); end for i=2:n for j=i:n z=0; for r=1:i-1 z=z+L(i,r)*U(r,j); end U(i,j)=A(i,j)-z; end if abs(U(i,i)) flag='failure' return; end for k=i+1:n m=0; for q=1:i-1 m=m+L(k,q)*U(q,i); end L(k,i)=(A(k,i)-m)/U(i,i); end end L U y=L\b;x=U\y detA=det(L*U) (2)列主元消去法: function x = gauss(A,b); A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]; b=[8;5.900001;5;1]; [n,n] = size(A); x = zeros(n,1); Aug = [A,b]; %增广矩阵 for k = 1:n-1 [piv,r] = max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r r = r + k - 1; % 列主元所在大矩阵的行 if r>k temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)==0, error(‘对角元出现0’), end % 把增广矩阵消元成为上三角 for p = k+1:n Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end end % 解上三角方程组 A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x(k)=x(k)/A(k,k); end detA=det(A) 2、运行结果 分解解线性方程组 1)LU L = 1.0e+006 * 0.0000 0 0 0 -0.0000 0.0000 0 0 0.0000 -2.5000 0.0000 0 0.0000 -2.4000 0.0000 0.0000 U = 1.0e+007 * 0.0000 -0.0000 0 0.0000 0 -0.0000 0.0000 0.0000 0 0 1.5000 0.5750 0 0 0 0.0000 x = -0.0000 -1.0000 1.0000 1.0000 detA = -762.0001 2)列主元消去法 detA = 762.0001 ans = 0.0000 -1.0000 1.0000 1.0000 3、问题结果 1)LU分解解线性方程组 1000,, ,,,3100,,L= ,,0.5250000010, ,,0.224000000.961,,, 10701,,, ,,0062.3,,U= ,,00150000055749998.499 ,,0005.079998908,, ,x=(,,,,,,，,,,,,,,，,,,,,,，,,,,,,) detA=-762.001 2)列主元消去法 ,x=(,,,,,,，,,,,,,,，,,,,,,，,,,,,,) detA=762.001 第六章 1、程序代码 (1)Jacobi迭代 Clear;clc; n = 6; %也可取n=8,10 H = hilb(n); b = H * ones(n, 1); e = 0.00001; for i = 1:n if H(i, i)==0 '对角元为零，不能求解' return end end x = zeros(n, 1); k = 0; kend = 10000; r = 1; while k<=kend &="" r="">e x0 = x; for i = 1:n s = 0; for j = 1:i - 1 s = s + H(i, j) * x0(j); end for j = i + 1:n s = s + H(i, j) * x0(j); end x(i) = b(i) / H(i, i) - s / H(i, i); end r = norm(x - x0, inf); k = k + 1; end if k>kend '迭代不收敛，失败' else '求解成功' x k end (2)SOR迭代 1)程序代码 function s = SOR(n, w); H = hilb(n); b = H*ones(n, 1); e = 0.00001; for i = 1:n if H(i,i)==0 ‘对角线为零，不能求解’ return end end x = zeros(n, 1); k = 0; kend = 10000; r = 1; while k<=kend &="" r="">e x0 = x; for i = 1:n s = 0; j = 1:i - 1 for s = s + H(i, j) * x(j); end for j = i + 1:n s = s + H(i, j) * x0(j); end x(i) = (1 - w) * x0(i) + w / H(i, i) * (b(i) - s); end r = norm(x - x0, inf); k = k + 1; end if k>kend '迭代不收敛，失败' else '求解成功' x end 2) 从命令窗口中分别输入: SOR(6,1) SOR(8,1) SOR(10,1) SOR(6,1.5) SOR(8,1.5) SOR(10，1.5) 2、运行结果 Jacobi迭代: ans = 迭代不收敛，失败 SOR迭代: 第七章 1、程序代码 (1)不动点迭代法 1)建立函数文件:g.m function f=g(x) f(1)=20/(x^2+2*x+10); 2)建立函数文件:bdd.m function [y, n] = bdd(x, eps) if nargin==1 eps=1.0e-8; elseif nargin<1> error return end x1 = g(x); n = 1; while (norm(x1-x)>=eps)&&(n<=10000)> x = x1; x1 = g(x); n = n + 1; end y = x; n 3)从命令窗口输入:bdd(0) (2)牛顿迭代 clear;clc; format long; m=8; %m为迭代次数，可分别令m=2,4,6,8,10 x=sym('x'); f=sym('x^3+2*x^2+10*x-20'); df=diff(f,x); FX=x-f/df; %牛顿迭代公式 Fx=inline(FX); disp('x='); x1=0.5; disp(x1); Eps=1E-8; k=0; while 1 x0=x1; k=k+1; x1=feval(Fx,x1); %将x1代入牛顿迭代公式替代x1 disp(x1); %在屏幕上显示x1 if k==m break; end end k,x1 2、运行结果 (1)不动点迭代法 >> bdd(0) n = 25 ans = 1.3688 (2) 牛顿迭代 x= 0 2 1.466666666666667 1.371512013805921 1.368810222633895 1.368808107822667 1.368808107821373 1.368808107821373 1.368808107821373 k = 8 x1 = 1.368808107821373 3、问题结果 (1)不动点迭代法 x=1.3688 n=25 收敛太慢 (2)牛顿迭代 初值取0 迭代次数k=8时， k x(k) 1 2 2 1.4666666 3 1.3715120 4 1.3688102 5 1.3688081 6 1.3688081 7 1.3688081 8 1.3688081 转载请注明出处范文大全网 » 数值计算(数值分析)试题及答范文四：数值分析计算实习题
范文五：数值分析计算实习题
