计算方法第七章

合集下载

地下结构的地层结构计算方法

开挖效应与释放荷载
模型建立要点
midas地层结构算例
第七章地层结构法的适用性
位移清零
模型建立要点
midas地层结构算例
初始地应力场
计算开挖边界等效结点力
删除开挖网格反向施加结点力
确定释放系数
第七章地层结构法的适用性
荷载分步释放与围岩特性
岩爆
模型建立要点
midas地层结构算例
高地应力
0.7m
E砼＝23Gpa A砼＝0.28m2 I砼＝0.00183m4
E钢＝210Gpa A钢＝39.578×10-4m2 I钢＝2500×10-8m4
E A ＝ E砼 A砼+ E钢A钢/S E I ＝ E砼 I砼+ E钢I钢/S
取E ＝ E砼
A ＝ A砼+ E钢A钢/(SE砼) =0.3316 I ＝ I砼+ E钢I钢/(SE砼) =0.002155
岩土材料
• 根据岩土性质和计算目的选择适合的本构模型。 • 定量分析时应注意材料参数的确定，必要时采用反分析。
结构材料
• 弹性或弹塑性 • 初期支护内的钢拱架与喷射砼一般视为整体计算
加固地层材料
• 直接模拟 • 不模拟，作为安全储备 • 提高地层材料参数
第七章地层结构法的适用性
边界条件
模型建立要点
576个四边形单元
35个梁单元
第七章地层结构法的适用性
模型建立要点
midas地层结构算例
地层与结构连接
公共节点，变形协调
. . . 1 node . A. B.
不同节点，相互独立
. . .. . 2 nodes . A. B.
摩擦接触，接触单元

计算方法第七章(特征值与特征向量)

( j p, q) i 1, 2, , n
最后，雅可比方法的计算步骤可以归纳为：（1）确定非对角绝对值最大元位置（p，q），并计算sin和 cos的值；（2）计算迭代矩阵的元素；
（3）计算特征向量；
（4）与计算精度进行比较，以决定第三节 QR 分解方法 3.1 QR 分解设 u 为n维实单位向量，称下面矩阵为Householder矩阵：
则
(2) (3) 1 a12 a13 (3) a 2 23 (3) Q2 A1 Q2Q1 A a33 (3) 0 a 3n
埃特金加速：可以证明：乘幂法线性收敛
mk 1 1
2 mk 1 1
2 1
[ zk 1 10 ] i [ zk 10 ] i
2 1
称为收敛率
由于
zk
线性收敛于 x1 ，于是可以对之进行埃特金加速，
( zk )i ( zk 2 )i ( zk 1 )i2 Wi ( zk )i 2( zk 1 )i ( zk 2 )i
, a
(k) pq
0
第 k 步迭代矩阵的元素为：
a a a
(k ) pj
a a
( k 1) pj
cos a
2
( k 1) qj
sin a
(k ) jp
(k ) k 1) ( k 1) k) aqj a (pj sin aqj cos a (jq ( j p, q ) (k ) pp ( k 1) pp
cos 2a a
( k 1) pp
(k 1) pq
sin cos a
( k 1) pq
(k 1) qq

电子教案第七章产品成本计算的基本方法——分批法

第七章产品成本计算的基本方法——分批法本章思维导图课题产品成本计算的分批法课时6课时（270min）教学目标知识技能目标：1.掌握分批法的特点、适用范围、计算程序和计算方法2.掌握简化分批法的应用条件、基本生产二级帐的作用以及在生产费用分配上的特点。

素质目标：培养学生的辩证思维和逻辑思维能力，培养学生在国家先进制造战略下将所学知识应用到实践问题分析中，做到学以致用。

教学重难点教学重点：分批法的概念及适用范围，简化分批法的概念和特点教学难点：分批法和简化分批法的成本计算程序教学方法讲授法、问答法、案例分析法、讨论法教学用具电脑、投影仪、多媒体课件、教材教学设计第1节课：考勤（2 min）→传授新知（28 min）→课堂讨论（15 min）第2节课：问题导入（5 min）→传授新知（20 min）→课堂练习（15 min）→课堂小结（3 min）→作业布置（2 min）教学过程主要教学内容及步骤设计意图第一节课问题导入【教师】提出问题：通过问题导入现代环境下怎样的企业适合用分批法法？与品种法有何区别？⏹【学生】思考、分组讨论并回答⏹【教师】通过学生的回答引入新课题：产品成本计算的分批法，引导学生思考分批法的概念，及其与品种法的区别，激发学生的学习兴趣传授新知⏹【教师】讲解分批法的概念、适用范围、特点、成本计算程序和简化分批法一、分批法的概念及适用范围分批法也称订单法，是以产品的批别或订单作为成本计算对象，归集和分配生产费用、计算产品成本的方法。

这种方法主要适用于小批、单件，管理上不要求分步骤计算成本的多步骤企业，如精密仪器、专用设备、重型机械和船舶制造企业，也可用于新产品的试制、机器设备修理、来料加工和辅助生产的工具模具制造企业等。

【教师】播放视频“分批法”（详见教材），帮助学生了解分批法【学生】观看、思考、理解二、分批法的特点（一）成本计算对象分批法是以产品的批别作为成本计算对象，这也是它区别于其他成本计算方法的最重要特征。

第七章结构位移计算

W=FP△ = FP△`cosa
第七章结构位移计算(Displacement)
2、静力实功在静外力FP1作用下，变形体在力的作用点沿力的方向发生位移△11 。静力实功为：式中的“1/2” ？ W=(1/2)FP1△11
静力概念：静力荷载加载到结构上是有一个过程的，在这个加载过程中，荷载从零增加到最后值，结构的内力和位移也达到最后值；在整个加载过程中，外力和内力始终保持静力平衡。
第七章结构位移计算(Displacement)
⒋ 本章在全课程中的地位想求静定结构的位移，必先求出静定结构的内力。因此本章可以说是对前面所学的各类静定结构的内力计算的复习与巩固。同时，位移计算又是下章即将开始学习的超静定结构的基础。因而,从全课程来看，本章是承上启下的一章，也是十分重要的内容。希望每个同学重视起来。
D
变形位移
ABDC ABD”C” 刚体位移
C
D
ABD’C’ 变形位移
位移状态
第七章结构位移计算(Displacement)
§7-2 变形体系的虚功原理
⒉ 着眼于位移：
dW ＝ dW dW dW ＝ dW
总总刚刚
变
微段平衡，由刚体虚功原理
dW刚 0
W总 W变
第七章结构位移计算(Displacement)
(a) 根据叠加原理，图(a)可分解为图(b)、(c)两种情况。 ※一个结构的两种状态。
(b)
(c)
第七章结构位移计算(Displacement) §7-2 变形体系的虚功原理一、刚体系的虚位移原理
刚体系处于平衡的充要条件是：在具有理想约束的
⒋ 用于动力计算和稳定性计算。

第七章土压力计算

粘性土主动土压力强度包括两部分
1. 土的自重引起的土压力zKp
2. 粘聚力c引起的侧压力2c√Kp 说明：侧压力是一种正压力，在计算 Ep 中应考虑
土压力合力
E p(1/2)h2K p2chK p
1.粘性土被动土压力强度不存在负侧压力区 2.合力大小为分布图形的面积，即梯形分布图形面积 3.合力作用点在梯形形心
被动朗肯状态
处于被动朗肯状态，σ3方向竖直，剪切
破坏面与竖直面夹角为45o＋/2
二、主动土压力
挡土墙在土压力作用下，产生离开
土体的位移，竖向应力保持不变，
水平应力逐渐减小，位移增大到
h
z
z(σ1)
-△a，墙后土体处于朗肯主动状态
时，墙后土体出现一组滑裂面，剪
45o＋/2
pa(σ3) 切破坏面与大主应力作用面夹角
z
pp
主动极限水平方向均匀伸展土体处于水平方向均匀压缩被动极限
平衡状态
弹性平衡
平衡状态
状态
主动朗肯状态
处于主动朗肯状态，σ1方向竖直，剪切
破坏面与竖直面夹角为45o-/2
被动朗肯状态
处于被动朗肯状态，σ3方向竖直，剪切
破坏面与竖直面夹角为45o＋/2
成层填土情况（以无粘性土为例）
A pa A
ppzK p2c Kp
pp zKp
h
h/3
Ep (1/2)h2Kp
hKp 1.无粘性土被动土压力强度与z成正比，沿墙高呈三角形分布 2.合力大小为分布图形的面积，即三角形面积 3.合力作用点在三角形形心，即作用在离墙底h/3处
h
hp
当c＞0, 粘性土
2c√Kp
hKp ＋2c√Kp

最优化计算方法

它的一般形式是：它的一般形式是：
min
f = c1x1 + c 2 x 2 + + c n x n a 11x1 + a 12 x 2 + + a1n x n <= b1 a x + a x + + a x <= b 21 1 22 2 2n n 2 a m1x1 + a m 2 x 2 + + a mn x n <= b m x i >= 0 (i = 1,2,, n )
在命令窗口输入：在命令窗口输入： x0=[0;0]; x=fminunc(‘fun703’,x0) 结果显示：结果显示： f =5.2979e-011 x =1.0673 0.1392 则非线性方程组的解为x1=1.0673,x2=0.1392。。则非线性方程组的解为
Matlab程序：程序：程序 ch703.m
第七章最优化计算方法
第一节线性方程组的应用
一、实验目的：实验目的：
1、了解线性规划问题及可行解、最优解的概念；、了解线性规划问题及可行解、 2、掌握Matlab软件关于求解线性规划的语句和方法。、掌握软件关于求解线性规划的语句和方法。软件关于求解线性规划的语句和方法
二、实验原理和方法：实验原理和方法：
迭代的基本思想和步骤大致可分为以下四步：迭代的基本思想和步骤大致可分为以下四步：
1) 2) 3）选取初始点x 0 , 并令k = 0; 得到x k 后，选取一个搜索方向P k , 使得沿着这个方向的目标函数f ( x)的值时下降的；由x k出发，沿P k 方向选取适当的步长λk , 使得 f ( x k + λk P k ) < f ( x k ) 由此得到下一个点x k +1 = x k + λk P k 4) 检验新得到的点x k +1是否满足精度要求的最优解。如果是，则结束运算；否则，令k = k + 1, 返回(2)继续迭代

工程水文学第七章流域产汇流计算

Rp.t （5） 0 0 0.07 0.13 0.51 1.10 0 11.15 26.73
Wt （6） 18.93 17.38 16.06 17.80 20.65 29.11 41.36 38.81 89.70
备注
（7）
Wm=120 mm Kw,t=1.0 Wt为一日开始时蓄水量（mm ）5月14 日开始时 18.93m m，是由前面计算得到的
NCWU
NCWU
月.日（1） 5.14 15 16 17 18 19 20 21 22
Pt （2） 0 0 3 4.2 10.3 15.1 0 63.2 56.8
Ew.t （3） 9.8 9.1 8.9 8.2 7.7 7.2 7.4 3.6 3.2
Et （4） 1.55 1.32 1.19 1.22 1.33 1.75 2.55 1.16 2.39
NCWU
两时段降雨： P1=49mm P2=81mm
(130mm)
降雨开始时： Pa=60mm 由 P1=49mm,查得R1=20.0mm。
（49mm)
由 P1 +P2=130mm, 查得R1+R2=80.0mm。则第二时段净雨为 R2=80-20=60mm
NCWU
二 P+Pa～R两变量相关图法
第7章流域产汇流计算
研究对象:
从定量上研究降雨形成径流的原理和计算方法，包括流域的产流计算和汇流计算。产流计算主要研究流域上降雨扣除植物截留、下渗、填洼等损失，转化为净雨过程的计算方法。汇流计算主要研究净雨沿地面和地下汇入河网，并经河网汇集形成流域出口断面径流过程的计算方法。
研究目的:
研究的流域产汇流计算是工程水文学中最基本的概念和方法之一，是以后学习由暴雨资料推求设计洪水，降雨径流预报等内容的基础。

数值计算方法第七章

第七章常微分方程初值问题数值解法7.1 引言科学技术中常常需要求解常微分方程的定解问题．这类问题最简单的形式，是本章将要着重考察的一阶方程的初值问题()()'00,y f x y y x y ⎧=⎪⎨=⎪⎩ ()()1.11.2我们只要，只要函数(),f x y 适当光滑—譬如关于y 满足利普希茨()Lipschitz 条件()(),,f x y f x y L y y -≤-. ()1.3 理论上就可以保证初值问题()1.1，()1.2的解()y y x =存在并且唯一.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法.所谓数值解法，就是寻求解()y x 在一系列离散节点121n n x x x x +<<<<<上的近似值121,,,,n n y y y y +.相邻两个节点的间距1n n n h x x +=-成为步长.今后如不特别说明，总是假定()1,2,i h h i ==为定数，这时节点为0,0,1,2,.n x x nh n =+= 初值问题()1.1，()1.2的数值解法有个基本特点，它们都采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进．描述这类算法，只要给出用已知信息,12,,n n n y y y --计算1n y +的递推公式.首先，要对方程()1.1离散化，建立求数值解的递推公式．一类是计算1n y +时只用到前一点的值n y ，称为单步法．另一类是用到1n y +前面k 点的值11,,,n n n k y y y --+称为k 步法．其次，要研究公式的局部截断误差和阶，数值解n y 与精确解()n y x 的误差估计及收敛性，还有递推公式的计算稳定性等问题．7.2 简单的数值方法与基本概念7.2.1 欧拉法与后退欧拉法我们知道，在xy 平面上，微分方程()1.1的解()y y x =称为它的积分曲线．积分曲线上一点(),x y 的切线斜率等于函数(),f x y 的值．如果按函数(),f x y 在xy 平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点方向相一致．基于上述几何解释，我们从初始点()000,p x y 出发，先依方向场在该点的方向推进到1x x =上一点1p ，然后再从1p 依方向场的方向推进到2x x =上一点2p ，循此前进做出一条折线012p p p （图91-）．一般地，设已做出该折线的顶点n p ，过(),n n n p x y 依方向场的方向再推进到()111,n n n p x y +++，显然两个定点n p ，1n p +的坐标有关系()11,,_n nn n n ny y f x y x x ++-=即()1,n n n ny y h f xy +=+ ()2.1 这就是著名的欧拉()Euler 公式．若初值0y 已知，则依公式()2.1可逐步算出()()10002111,,,,y y hf x y y y hf x y =+=+例1 求解初值问题()'201x y y y y ⎧=-⎪⎨⎪=⎩()01x <<， ()2.2解为便于进行比较，本章将用多种数值方法求解上述初值问题．这里先用欧拉方法，欧拉公式的具体形式为12.n n n n n x y y h y y +⎛⎫=+- ⎪⎝⎭取步长0.1,h =计算结果见表71-表71-计算结果对比n xn y()n y xn xn y()n y x0.10.2 0.30.4 0.51.1000 1.1918 1.2774 1.3582 1.43511.0954 1.1832 1.2649 1.3416 1.41420.6 0.7 0.80.9 1.01.5090 1.5803 1.6498 1.7178 1.78481.4832 1.5492 1.6125 1.6733 1.7321初值问题()2.2有解12y x =+，按这个解析式子算出的准确值()n y x 同近似值n y 一起列在表71-中，两者相比较可以看出欧拉方法的精度很差．还可以通过几何直观来考察欧拉方法的精度．假设()n n y y x =，即顶点n p 落在积分曲线()y y x =上，那么，按欧拉方法做出的切线n p 1n p +便是()y y x =过点n p 的切线．从图形上看，这样定出顶点1n p +显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的．为了分析计算公式的精度，通常可采用泰勒展开将()1n y x +在n x 处展开，则有()()1n n y x y x h +=+ ()()()2''',2n n n h y x yx h y ξ=++ ()1,.n n n x x ξ+∈在()n n y y x =的前提下，()()()()',,.n n n n n f x y f x y x y x ==于是可得欧拉法()2.1的公式误差()()()22''''1122n n n n h h y x y y y x ξ++-=≈， ()2.3 称为此方法的局部截断误差．如果对方程()1.1从n x 到1n x +积分，得 ()()()()11,n n n nx y x y x f t y t dt x ++=+⎰． ()2.4 右端积分用左矩形公式()(),n n hf x y x 近似，再以n y 代替()n y x ，1n y +代替()1n y x +也得到()2.1，局部截断误差也是()2.3．如果在()2.4中右端积分用右矩形公式()()11,n n hf x y x ++近似，则得另一个式 ()()111,n n n n y y h f x y x +++=+ ()2.5称为后退得欧拉法．后退的欧拉公式与欧拉公式有着本质的区别，后者是关于1n y +的一个直接的计算公式，这类公式称作是显式的；然而公式()2.5的右端含有未知的1n y +，它实际上是关于1n y +的一个函数方程，这类公式称作是隐式的．显式与隐式两类方法各有特点．考虑到数值稳定性等其他因素，人们有时需要选用隐式方法，但使用显式算法远比隐式方便．隐式方程()2.5通常用迭代法求解，而迭代过程的实质是逐步显示化．设用欧拉公式()()01,n n n n y y hf x y +=+给出迭代初值()01n y +，用它代入()2.5式的右端，使之转化为显示，直接计算得()()()10111,n n n n y y hf x y +++=+，然后再用11n y +代入()2.5式，又有()()()21111,n n n n y y hf x y +++=+．如此反复进行，得()()()1111,k k n n n n y y hf x y ++++=+ ()0,1,k =． ()2.6由于(),f x y 对y 满足利普希茨()Lipschitz 条件()1.3．由()2.6减()2.5得()()()()1111111,,k k n n n n n n y y h f x y f x y +++++++-=-()11.k n n hL y y ++≤-由此可知，只要1hL <迭代法()2.6就收敛到解1n y +．关于后退欧拉法的公式误差，从积分公式看到它与欧拉法是相似的．7.2.2 梯形方法为得到比欧拉法精度高的计算公式，在等式()2.4右端积分中若用梯形求积公式近似，若用n y 代替()n y x ，1n y +代替()1n y x +，则得()()111,,2n n n n n n hy y f x y f x y +++=++⎡⎤⎣⎦ ()2.7 称为梯形方法．梯形方法是隐式单步法，可用迭代法求解．同后退的欧拉法一样，仍用欧拉方法提供迭代初值，则梯形法的迭代公式为()()()()()()()011111,;,,;20,1,2,.nn n n k k n n n n n n y y hf x y h y y f x y f x y k +++++⎧=+⎪⎪⎡⎤=++⎨⎣⎦⎪⎪=⎩()2.8 为了分析迭代过程的收敛性，将()2.7式与()2.8式相减，得()()()()1111111,,2k k n n n n n n h y y f x y f x y +++++++⎡⎤-=-⎣⎦，于是有()()11111,2k k n n n n hL y y y y +++++-≤- 式中L 为(),f x y 关于y 的利普希茨常数．如果选取h 充分小，使得12hL≤，则当k →∞时有()1k n y +1n y +→，这说明迭代过程()2.8是收敛的．7.2.3 单步法的局部截断误差与阶初值问题()()1.1,1.2的单步法可用一般形式表示为()11,,,,n n n n n y y h x y y h ϕ++=+ ()2.9其中多元函数ϕ与(),f x y 有关，当ϕ含有1n y +时，方法是隐式的，若不含1n y +则为显式方法，所以显式单步法可表示为 ()1,,,n n nn y y h x y hϕ+=+ ()2.10 (),,x y h ϕ称为增量函数，例如对欧拉法()2.1有(),,x y h ϕ(),f x y =．它的局部截断误差已由()2.3给出，对一般显式单步法则可如下定义．定义 1 设()y x 是初值问题()()1.1,1.2的准确解，称()()()()11,,n n n n n T y x y x h x y x h ϕ++=-- ()2.11为显式单步法()2.10的局部截断误差．1n T +之所以称为局部的，是假设在n x 前各步没有误差．当n y ()n y x =时，计算一步，则有()()()111,,n n n n n n y x y y x y h x y h ϕ+++-=-+⎡⎤⎣⎦()()()()11,,n n n n n y x y x h x y x h T ϕ++=--=．所以，局部截断误差可理解为用方法()2.10计算一步的误差，也即公式()2.10中用准确解()y x 代替数值解产生的公式误差．根据定义，显然欧拉法的局部截断误差()()()()11,n n n n n T y x y x h x y x ϕ++=--()()()'1n n n y x y x h y x +=--()()2''3,2n h y x h =+O 即为()2.3的结果．这里()2''2n h y x 称为局部截断误差主项．显然()21n T h +=O ．一般情形的定义如下，定义 2 设()y x 是初值问题()()1.1,1.2的准确解，若存在最大整数p 使显式单步法()2.10的局部截断误差满足()()()()11,,p n T y x h y x h x y h h ϕ++=+--=O ， ()2.12则称方法()2.10具有p 阶精度．若将()2.12展开式写成()()()121,p p n n n T x y x h h ψ+++=+O ，则()()1,p n n x y x h ψ+称为局部截断误差主项．以上定义对隐式单步法()2.9也是适用的．例如，对后退欧拉法()2.5其局部截断误差为()()()()1111,n n n n n T y x y x hf x y x ++++=--()()()2'''32n n h hy x y x h =++O()()()'''2n n h y x hy x h ⎡⎤-++O ⎣⎦()()2''3.2n h y x h =-+O这里1p =，是1阶方法，局部截断误差主项为()2''2n h y x -．同样对梯形法()2.7有()()()()''1112n n n n n h T y x y x y x y x +++⎡⎤=--+⎣⎦ ()()()23''''''23!n n n h h hy x y x y x =++()()()()()2'''''''422n n n n h h y x y x hy x y x h ⎡⎤-++++O ⎢⎥⎣⎦()()3'''4.12n h y x h =-+O所以梯形方法()2.7是二阶的，其局部误差主项为()3'''12n h y x -．7.2.4 改进的欧拉公式我们看到，梯形方法虽然提高了精度，但其算法复杂，在应用迭代公式()2.9进行实际计算时，每迭代一次，都要重新计算函数(),f x y 的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测．为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法．具体地说，我们先用欧拉公式求得一个初步的近似值1n y +，称之为预测值，预测值1n y +的精度可能很差，再用梯形公式()2.7将它校正一次，即按()2.8式迭代一次得1n y +，这个结果称校正值，而这样建立的预测-校正系统通常称为改进的欧拉公式：预测 1n y +(),,n n n y hf x y =+校正 1n y +()()11,,.2n n n nn hy f x y fx y ++⎡⎤=++⎣⎦()2.13或表为下列平均化形式()()()11,,,,1.2p n n n c n n p n p c y y hf x y y y hf x y y y y ++⎧=+⎪⎪⎪=+⎨⎪⎪=+⎪⎩ 例2 用改进的欧拉方法求解初值问题()2.2．解改进的欧拉公式为()112,2,1.2n p nn n n c n p p n p c x y y hf y y x y y hf y y y y y ++⎧⎛⎫=+-⎪ ⎪⎝⎭⎪⎪⎛⎫⎪=+- ⎪⎨ ⎪⎪⎝⎭⎪⎪=+⎪⎩仍取0.1h =，计算结果见表92-．同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度．表72- 计算结果对比n xn y()n y xn xn y()n y x0.10.2 0.30.4 0.51.0959 1.1841 1.2662 1.3424 1.41641.0954 1.1832 1.2649 1.3416 1.41420.6 0.7 0.80.9 1.01.4860 1.5525 1.6153 1.6782 1.73791.4832 1.5492 1.6165 1.6733 1.73217.3 龙格-库塔方法7.3.1 显式龙格-库塔法的一般形式上节给出了显式单步法的表达式()2.10，其局部截断误差为()2.12，对欧拉法()21n T h +=O ，即方法为1p =阶，若用改进欧拉法()2.13，它可表为()()()1,,,2n n n n n n n n hy y f x y f x h y hf x y +⎡⎤=++++⎣⎦． ()3.1 此时增量函数()()()()1,,,,,2n n n n n n n n x y h f x y f x h y hf x y ϕ⎡⎤=+++⎣⎦． ()3.2 它比欧拉法的(),,n n x y h ϕ= (),n n f x y ，增加了计算一个右函数f 的值，可望2p =．若要使得到的公式阶数p 更大，ϕ就必须包含更多的f 值．实际上从方程()1.1等价的积分形式()2.4，即()()()()11,n nx n n x y x y x f x y x dx ++-=⎰， ()3.3若要使公式阶数提高，就必须使右端积分得数值求积公式精度提高，它必然要增加求积节点，为此可将()3.3的右端用积分公式表示为()()()()11,,.n nrx i n i n i x i f x y x dx c f x h y x h λλ+=≈++∑⎰一般来说，点数r 越多，精度越高，上式右端相当于增量函数(),,x y h ϕ，为得到便于计算的显示方法，可类似于改进欧拉法()()3.1,3.2，将公式表示为 ()1,,,n n n n y y h x y h ϕ+=+ ()3.4其中()1,,rn n i i i x y h c K ϕ==∑， ()3.5(),,i n n K f x y =11,i i n i n ij j j K f x h y K λμ-=⎛⎫=++ ⎪⎝⎭∑ 2,,,i r = 这里i c ，i λ，ij μ均为常数．()3.4和()3.5称为r 级显式龙格-库塔()Runge Kutta -法，简称R K -方法．当1r =，()(),,,n n n n x y h f x y ϕ=时，就是欧拉法，此时方法的阶为1p =．当2r =时，要改进欧拉法()3.1，()3.2就是其中的一种，下面将证明阶2p =．要使公式()3.4，()3.5具有更高的阶p ，就要增加点数r ．下面我们只就2r =推导R K -方法．并给出3,4r =时的常用公式，其推导方法与2r =时类似，只是计算较复杂．7.3.2 二阶显式R K -方法对2r =的R K -方法，由()3.4，()3.5可得到如下的计算公式()()()11122122211,,,,.n n n n n n y y h c K c K K f x y K f x h y hK λμ+=++⎧⎪=⎨⎪=++⎩ ()3.6这里1c ，2c ，2λ，21μ均为待定常数，我们希望适当选取这些系数，使公式阶数p 尽量高．根据局部截断误差定义，()3.6的局部截断误差为()()11n n n T y x y x ++=-=()()12221,,,n nn n n h cf x y c f xh y h f λμ-+++⎡⎤⎣⎦()3.7 这里()(),,.n n n n n y y x f f x y ==为得到1n T +的阶p ，要将上式各项在(),n n x y 处做泰勒展开，由于(),f x y 是二元函数，估要用到二元泰勒展开，各项展开式为()()23''''''41,23!n n nn n h h y x y hy y y h +=++++O其中()()()()()()()()()()()'''''''''''2''''',,,,,,,2,,,,,;n n n n n n n x n n y n n n n xx n n n xy n n n yy n n y n n x n n n y n n y f x y f d y f x y x f x y f x y f dx y f x y f f x y f f x y f x y f x y f f x y ⎧==⎪⎪==+⋅⎨⎪⎪⎡⎤=++++⎣⎦⎩()3.8()221,n n n f x h y hf λμ++()()()''2221,,.n x n n y n n n f f x y h f x y hf h λμ=+++O 将以上结果代入()3.7则有()()2''1,,2n n x n n y n n n h T hf f x y f x y f +⎡⎤=++⎣⎦ ()()()()''312221,,n nx n n y n n n h c f c f f x y h f x y f h O h λμ⎡⎤-++++⎣⎦ ()()'2122211,2n x n n c c f h c f x y h λ⎛⎫=--+- ⎪⎝⎭()()'232211,2y n n n c f x y f h O h μ⎛⎫+-+ ⎪⎝⎭．要使公式()3.6具有2p =阶，必须使12222211110,0,0.22c c c c λμ--=-=-= ()3.9即222211211,, 1.22c c c c λμ==+=()3.9的解是不唯一的．可令20c a =≠，则得122111,.2c a aλμ=-==这样得到的公式称为二阶R-K 方法，如取1/2a =，则121/2c c ==，2211λμ==．这就是改进欧拉法()3.1．若取1a =，则211,0,c c ==2211/2λμ==．得计算公式()12121,,,,.22n n n n n n y y hK K f x y h h K f x y K +⎧⎪=+⎪⎪=⎨⎪⎛⎫⎪=++ ⎪⎪⎝⎭⎩()3.10称为中点公式，相当于数值积分公式的中矩形公式．()3.10也可表示为()1,,.22n n n n n n h h y y hf x y f x y +⎛⎫=+++ ⎪⎝⎭对2r =的R-K 公式()3.6能否使局部误差提高到()4O h ？为此需把2K 多展开一项，从()3.8的'''ny 看到展开式中()2'''yxyf f ff +的项是不能通过选择参数消掉的，实际上要使3h 的项为零，需增加3个方程，要确定4个参数122,,c c λ及21μ，这是不可能的．故2r =的显示方法的阶只能是2p =，而不能得到三阶公式．7.3.3 三阶与四阶显式R K -方法要得到三阶显式R K -方法，必须3r =．此时()()3.4,3.5的公式表示为()()()()111223312221133311322,,,,,,.n n n n n n n y y h c K c K c K K f x y K f x h y hK K f x h hK hK λμλμμ+=+++⎛= =++ =++⎝ ()3.11其中123,,c c c 及22133132,,,,λμλμμ均为待定参数，公式()3.11的局部截断误差为()()[]11112233.n n n T y x y x h c K c K c K ++=--++只要将23,K K 按二元函数泰勒展开，使()41n T O h +=，可得待定参数满足方程12322133132223322223332321,,,1,21,31.6c c c c c c c c λμλμμλλλλλμ++=⎧⎪=⎪⎪=+⎪⎪+=⎨⎪⎪+=⎪⎪⎪=⎩()3.12这是8个未知数6个方程的方程组，解也不是唯一的．可以得到很多公式．满足条件的公式统称为三阶R K -公式．下面只给出其中一个常见的公式．()()()11231213124,6,,,,22,2.n n n n n n n n h y y K K K K f x y h h K f x y K K f x h y hK hK +⎧=+++⎪⎪=⎪⎨⎛⎫⎪=++ ⎪⎪⎝⎭⎪=+-+⎩ 此公式称为库塔三阶方法．继续上述过程，经过较复杂的数学演算，可以导出各种四阶龙格-库塔公式，下列经典公式是其中常用的一个：()()()11234121324322,6,,,,.22,,22,.n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪⎪=⎪⎪⎛⎫=++⎨ ⎪⎝⎭⎪⎪⎛⎫=++⎪ ⎪⎝⎭⎪⎪=++⎩ ()3.13四阶龙格-库塔方法的每一步需要四次函数值f ，可以证明其截断误差为()5O h ．不过证明极其繁琐，这里从略．例3 设取步长0.2h =，从0x =直到1x =用四阶龙格-库塔方法求解初值问题()2.2．解这里，经典的四阶龙格-库塔公式()3,13具有形式右面列出计算结果n y ，表中()n y x 仍表示准确解．比较例3和例2的计算结果，显然以龙格-库表73- 计算结果塔方法的精度为高．要注意，虽然四阶龙格-库塔方法的计算量（每一步要4次计算函数f ）比改进的欧拉方法（它是一种二阶龙格-库塔方法，每一步只要2次计算函数f ）打一倍，但由于这里放大了步长，表和表所耗费的计算量几乎相同．这个例子又一次显示了选择算法的重要意义．然后值得指出的是，龙格-库塔方法的推导基于泰勒展开方法，因而它要求所求的解具有较好的光滑性质．反之，如果解得光滑性差，那么，使用四阶龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法．实际计算时，我们应当针对问题的具体特点选择合适的算法．7.3.4 变步长的龙格-库塔方法单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求n x n y ()n y x0.2 0.4 0.6 0.8 1.01.1832 1.3417 1.4833 1.6125 1.73211.18321.3417 1.4833 1.6125 1.7321()()11234121132243322,62,2,222,222.n n n n n n n n n n n n n n h y y K K K K x K y y x h h K y K h y K x h h K y K h y K x h K y hK y hK +⎧=++++⎪⎪⎪=-⎪⎪+⎪=+-⎪+⎨⎪+⎪=+-⎪+⎪⎪⎪+=+-⎪+⎩解范围内所要完成的步数就增加了．步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累．因此同积分的数值计算一样，微分发出的数值解法也有个选择步长的问题．在选择步长时，需要考虑两个问题：1 怎样衡量和检验计算结果的精度？2 如何依据所获得的精度处理步长？我们考察经典的四阶龙格-库塔公式()3,13．从节点n x 出发，先以h 为步长求出一个近似值，记为()1h n y +，由于公式的局部截断误差为()5O h ，故有 ()()511h n n y x y ch ++-≈， ()3.14 然后将步长折半，即取2h 为步长从n x 跨两步到1n x +，再求得一个近似值21h n y ⎛⎫⎪⎝⎭+，每跨一步的截断误差是52h c ⎛⎫⎪⎝⎭，因此有()521122h n n h y x yc ⎛⎫ ⎪⎝⎭++⎛⎫-≈ ⎪⎝⎭， ()3.15比较()3.14式和()3.15式我们看到，步长折半后，误差大约减少到116，即有 ()()()211111.16h n n h n n y x yy x y⎛⎫ ⎪⎝⎭++++-≈- 由此易得下列事后估计式()()2211111.15h h h n n n n y x yy y ⎛⎫⎛⎫⎪ ⎪⎝⎭⎝⎭++++⎡⎤-≈-⎢⎥⎢⎥⎣⎦这样，我们可以通过检查步长，折半前后两次计算结果的偏差()211h h n n yy ⎛⎫ ⎪⎝⎭++∆=-来判定所选的步长是否合适，具体地说，将区分以下两种情况处理：⒈ 对于给定的精度ε，如果ε∆>，我们反复将步长折半进行计算，直至ε∆<为止，这样取最终得到的21h n y⎛⎫ ⎪⎝⎭+作为结果；⒉ 如果ε∆<，我们将反复将步长加倍，直到ε∆>位止，这时再将步长折半一次，就得到所要的结果．这种通过加倍或折半步长的方法称为变步长方法．表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的．7.4 单步法的收敛性与稳定性7.4.1 收敛性与相容性数值解法的基本思想是，通过某种离散化手段将微分方程()1.1转化为差分方程，如单步法()2.10，即()1,,.n n nn y y h x y h ϕ+=+ ()4.1它在n x 处的解为n y ，而初值问题()()1.1,1.2在n x 处的精确解为()n y x ，记()n n n e y x y =-称为整体截断误差．收敛性就是讨论当n x x =固定且00n x x h n-=→时0n e →的问题．定义3 若一种数值方法（如单步法()4.1）对于固定的0n x x nh =+，当0h →时有()n n y y x →，其中()y x 是()()1.1,1.2的准确解，则称该方法是收敛的．显然数值方法收敛是指()0n n n e y x y =-→，对单步法()4.1有下述收敛定理：定理1 假设单步法()4.1具有p 阶精度，且增量函数(),,x y h ϕ关于y 满足利普希茨条件(),,,,x y h x y h L y y ϕϕϕ--⎛⎫-≤- ⎪⎝⎭， ()4.2又设是准确的，即()00y y x =，则其整体截断误差()().p n n y x y O h -= ()4.3证明设以1n y -+表示取()n n y y x =公式()4.1求得的结果，即()()()1,,n n n n y y x h x y x h ϕ-+=+， ()4.4则()11n n y x y -++-为局部截断误差，由于所给方法具有p 阶精度，按定义2，存在定数C ，使()111p n n y x y Ch -+++-≤．又由式()4.4与()4.1，得()11n n n n y y y x y -++-≤- ()()(),,,,.n n n nh x yx h x y h ϕϕ+- 利用假设条件()4.2，有()()111,n n n n y y hL y x y ϕ-++-≤+-从而有()()111111n n n n n n y x y y y y x y --++++++-≤-+-()()11.p n n hL y x y Ch ϕ+≤+-+即对整体截断误差()n n n e y x y =-成立下列递推关系式()111p n n e hL e Ch ϕ++≤++， ()4.5据此不等式反复递推，可得()()0111.p nn n Ch e hL e hL L ϕϕϕ⎡⎤≤+++-⎢⎥⎣⎦ ()4.6 再注意到当0n x x nh T -=≤时()()1nnhL TL hL eeϕϕϕ+≤≤，最终得下列估计式()01.p TL TL n Ch e e ee L ϕϕϕ≤+- ()4.7由此可以断定，如果初值是准确的，即00e =，则()4.3式成立，定理证毕．依据这一定理，判断单步法()4.1的收敛性，归结为验证增量函数能否满足利普希茨条件．对于欧拉方法，由于其增量函数ϕ就是(),f x y ，故当(),f x y 关于y 满足利普希茨条件时它是收敛的．再考察改进的欧拉方法，其增量函数已由()3.2式给出，这时有()()1,,,,[,(,)2x y h x y h f x y f x y ϕϕ--⎛⎫-≤- ⎪⎝⎭((,(,))(,(,))]f x h y h f x y f x h y h f x y --+++-++．假设(),f x y 关于y 满足利普希茨条件，记利普希茨常数L ，则由上式推得(),,,,1.2h x y h x y h L L y y ϕϕ--⎛⎫⎛⎫-≤+- ⎪ ⎪⎝⎭⎝⎭设限定0h h ≤（0h 为定数）上式表明ϕ关于y 的利普希茨常数01,2h L L L ϕ⎛⎫=+ ⎪⎝⎭因此改进的欧拉方法也是收敛的．类似地，不难验证其他龙格-库塔方法的收敛性．定理1表明时单步法收敛，并且当是初值问题的解，具有阶精度时，则有展开式1()()(,(),)n T y x h y x h x y x h ϕ+=+--'''2()()2y x y x h h =++'[(,(),0)(,(),0)]xh x y x x y x h ϕϕ-++ '2[()(,(),)]().h y x x y x h O h ϕ=-+所以1p ≥的充要条件是'()(,(),0)0y x x y x ϕ-=，而'()(,())y x f x y x =，于是可给出如下定义：定义4 若单步法()4.1的增量函数ϕ满足(,,0)(,)x y f x y ϕ=，则称单步法()4.1与初值问题()()1.1,1.2相容．以上讨论表明p 阶方法()4.1当1p ≥时与()()1.1,1.2相容，反之相容方法至少是1阶的．于是由定理1可知方法()4.1收敛的充分必要条件是此方法是相容的．7.4.2 绝对稳定性与绝对稳定前面关于收敛性的讨论有个前提，必须假定数值方法本身的计算是准确的．实际情形并不是这样，差分方程的求解还会有计算误差，譬如由于数字舍入而引起的小扰动．这类小扰动在传播过程中会不会恶性增长，以至于“淹没”了差分方程的“真解”呢？这就是差分方法的稳定性问题．在实际计算时，我们希望某一步产生的扰动值，在后面的计算中能够被控制，甚至是逐步衰减的．定义5 若一种数值方法在节点值n y 上大小为δ的扰动，于以后各节点值()m y m n >上产生的误差均不超过δ，则称该方法是稳定的．下面先以欧拉法为例考察计算稳定性．例4 考察初值问题()'100,0 1.y y y ⎧=-⎪⎨=⎪⎩ 其准确解()100x y x e -=是一个按指数曲线衰减得很快的函数，用欧拉法解方程'100y y =- 得()11100n n y h y +=-．若取0.025h =，则欧拉公式的具体形式为1 1.5,n n y y +=-计算结果列于表94-的第2列．我们看到，欧拉方法的解n y 在准确值()n y x 的上下波动，计算过程明显地不稳定．但若取10.005,0.5n n h y y +==则计算过程稳定．在考察后退的欧拉方法，取0.025h =时计算公式为113.5n n y y +=．计算结果列于表74-的第3列，这时计算过程是稳定的．表74- 计算结果对比节点欧拉方法后退欧拉方法 0.025 0.050 0.075 0.1001.5-2.253.375- 5.06250.2857 0.0816 0.0233 0.0067例题表明稳定性不但与方法有关，也与步长h 的大小有关，当然也与方程中的(),f x y 有关．为了只考察数值方法本身．通常只检验将数值方法用于解模型方程的稳定性，模型方程为'y y λ=， ()4.8其中λ为复数，这个方程分析较简单．对一般方程可以通过局部线性化化方程为这种形式，例如在(,)x y --的邻域，可展开为()',y f x y =''(,)(,)()(,)(),x y f x y f x y x x f x y y y --------=+-+-+略去高阶项，再做变换即可得到'u u λ=的形式．对于m 个方程的方程组，可线性化为'y Ay =，这里A 为m m ⨯的雅可比矩阵ii fy ⎛⎫∂ ⎪∂⎝⎭．若A 有m 个特征值1,,m λλ，其中i λ可能是复数．所以，为了使模型方程结果能推广到方程组，方程()4.8中λ为复数．为保证微分方程本身的稳定性，还应假定()Re 0λ<．下面先研究欧拉方法的稳定性．模型方程'y y λ=的欧拉公式为()11.n n y h y λ+=+ ()4.9设在节点值n y 上有一扰动值n ε，它的传播使节点值1n y +产生大小为1n ε+的扰动值，假设用*n n n y y ε=+按欧拉公式得出*111n n n y y ε+++=+的计算过程不再有新的误差，则扰动值满足()11.n n h ελε+=+可见扰动值满足原来的差分方程()4.9．这样，如果差分方程的解是不增长的，即有1n n y y +≤，则它就是稳定的．这一论断对于下面将要研究的其他方法同样适用．显然，为要保证差分方程()4.9的解是不增长的，只要选取h 充分小，使1 1.h λ+≤ ()4.10在h μλ=的复平面上，这是以()1,0-为圆心，1为半径的单位圆域．称为欧拉法的绝对稳定域，一般情形可由下面定义．定义6 单步法()4.1用于解模型方程()4.8，若得到的解()1n n y E h y λ+=，满足()1E h λ<，则称方法()4.1是绝对稳定的．在h μλ=的平面上，使()1E h λ<的变量围成的区域，称为绝对稳定域，它与实轴的交称为绝对稳定区间．对欧拉法()1E h h λλ=+，其绝对稳定域已由()4.10给出，绝对稳定区间为20h λ-<<，在例5中100,21000h λ=--<-<，即02/1000.02h <<=为绝对稳定区间，例4中取0.025h =故它是不稳定的，当取0.005h =时它是稳定的．对二阶R K -方法，解模型方程()4.8可得到()211,2n n h y h y λλ+⎡⎤=++⎢⎥⎢⎥⎣⎦故()()21.2h E h h λλλ=++绝对稳定域由()2112h h λλ++<得到，于是可得到绝对稳定区间为20h λ-<<，即02/h λ<<-．类似可得到三阶及四阶的R K -方法的()E h λ分别为()()()2312!3!h h E h h λλλλ=+++，()()()()2341.2!3!4!h h h E h h λλλλλ=++++由()1E h λ<可得到相应的绝对稳定域．当λ为实数时则得绝对稳定区间．它们分别为三阶R K -显示方法： 2.510,h λ-<<即0 2.51/.h λ<<- 四阶R K -显示方法： 2.780,h λ-<<即0 2.78/h λ<<-．从以上讨论可知显示的R K -方法的绝对稳定域均为有限域，都对步长h 有限制．如果h 不在所给的绝对稳定区间内，方法就不稳定．例5 '20y y =- ()01x ≤≤，()01y =分别取0.1h =及0.2h =用经典的R K -四阶方法()3.13计算．解本例20λ=-，h λ分别为2-及4-，前者在绝对稳定区间内，后者则不在，用四阶R K -方法计算其误差见下表：n x 0.20.40.60.81.00.1h = 0.2h =10.9310-⨯ 4.98 10.1210-⨯ 25.0 20.1410-⨯ 125.0 30.1510-⨯ 625.0 40.1710-⨯3125.0从以上结果看到，如果步长h 不满足绝对稳定条件，误差增长很快．对隐式单步法，可以同样讨论方法的稳定性，例如对后退欧拉法，用它解模型方程可得11,1n n y y h λ+=- 故 ()11E h h λλ=-．由()111E h h λλ=<- 可得绝对稳定域为11h λ->，它是以()1,0为圆心，1为半径的单位圆外部，故绝对稳定区间为0h λ-∞<<．当0λ<时，则0h <<∞，即对任何步长均为稳定的．对隐式梯形法，它用于解模型方程()4.8得112,12n n h y y h λλ++=-故 ()1/2.1/2h E h h λλλ+=- 对()Re 0λ<有()12112h E h h λλλ+=<-，故绝对稳定域为h μλ=的左半平面，绝对稳定区间为0h λ-∞<<，即0h <<∞时梯形法均是稳定的．7.5 线性多步法在逐步推进的求解过程中，计算1n y +之前事实上已经求出了一系列的近似值01,,,n y y y ，如果充分利用前面多步的信息来预测，则可以期望会获得较高的精度．这就是构造所谓线性多步法的基本思想．构造多步法的主要途径是基于数值积分方法和基于泰勒展开方法，前者可直接由方程()1.1两端积分后利用插值求积公式得到．本节主要介绍基于泰勒展开的构造方法．7.5.1 线性多步法的一般公式如果计算n k y +时，除用1n k y +-的值，还用到()0,1,,2n i y i k +=-的值，则称此方法为线性多步法．一般的线性多步法公式可表示为10,k kn k i n i i n i i i y y h f αβ-+++===+∑∑ ()5.1其中n i y +为()n i y x +的近似，()0,,,,n i n i n i n i i i f f x y x x ih αβ++++==+为常数，0α及0β不全为零，则称()5.1为线性k 步法，计算时需先给出前面k 个近似值011,,,k y y y -，再由()5.1逐次求出1,,k k y y +．如果0k β=，称()5.1为显示k 步法，这时n k y +可直接由()5.1算出;如果0k β≠，则()5.1称为隐式k 步法，求解时与梯形法()2.7相同，要用迭代方法可算出n k y +．()5.1中系数i α及i β可根据方法的局部截断误差及阶确定，其定义为：定义7 设()y x 是初值问题()()1.1,1.2的准确解，线性多步法()5.1在n k x +上的局部截断误差为();n k n T L y x h +=⎡⎤⎣⎦()()()1'0.k kn k i n i i n i i i y x y x h y x αβ-+++===--∑∑ ()5.2若()1p n k T O h ++=，则称方法()5.1是p 阶的，1p ≥则称方法()5.1与方程()1.1是相容的．由定义7，对n k T +在n x 处做泰勒展开，由于()()()()()2'''2!n n n n ih y x ih y x ihy x y x +=++()()3'''3!n ih y x ++，()()()()()2'''''''2!n n n n ih y x ih y x ihy x y x +=+++代入()5.2得()()()'2''012n k n n n T c y x c hy x c h y x +=++ ()()p p p n c h y x +++， ()5.3其中()0011,k c αα-=-++()()11210121,k k c k k αααβββ-=-+++--+++⎡⎤⎣⎦()()121121!qqq p k c k k q ααα-⎡⎤=-+++-⎣⎦()1112121!q q k k q βββ--⎡⎤-+++⎣⎦-2,3,.q = ()5.4若在公式()5.1中选择系数i α及i β，使它满足010,p c c c ==== 10p c +≠由定义可知此时所构造的多步法是p 阶的，且()()()1121.p p p n k p n T c h y x O h +++++=+ ()5.5称右端第一项为局部截断误差主项，1p c +称为误差常数．根据相容性定义，1p ≥，即010c c ==，由1,n n n y y hf +=+()5.4得0111101,.k k k i ii i i k ααααβ--==+++=⎧⎪⎨+=⎪⎩∑∑ ()5.6故方法与微分方程相容的充分必要条件是成立．显然，当1k =时，若10β=，则由()5.6可求得001, 1.αβ==此时公式()5.1为1,n n n y y hf +=+即为欧拉法．从()5.4可求得21/20c =≠，故方法为1阶精度，且局部截断误差为()()2''3112n n T h y x O h +=+，这和第2节给出的定义及结果是一致的．对1k =，若10β≠，此时方法为隐式公式，为了确定系数001,,αββ，可由0120c c c ===解得0011,1/2αββ===．于是得到公式()112n n n n hy y f f ++=++，即为梯形法．由()5.4可求得31/12c =-，故2p =，所以梯形法是二阶方法，其局部截断误差主项是()3'''/12n h y x -，这与第2节中的讨论也是一致的．对2k ≥的多步法公式都可利用()5.4确定系数,i i αβ，并由()5.5给出局部截断误差，下面只就若干常用的多步法导出具体公式．7.5.2 米尔尼方法与辛普森方法考虑另一个4k =的显示公式()13322110,n n n n n ny y h f f f f ββββ++++=++++ 其中0123,,,ββββ为待定常数，可根据使公式的阶尽可能高这一条件来确定其数值．由()5.4可知00c =，再令12340c c c c ====得到 ()()()01231231231234,22316,34964,4827256.βββββββββββββ+++=⎧⎪++=⎪⎨++=⎪⎪++=⎩ 解此方程组得3210848,,,0333ββββ==-== 于是得到四步显示公式()4321422,3n n n n n hy y f f f ++++=+-+ ()5.11 称为米尔尼（Milne ）方法.由于514/45c =，故方法为4阶，其局部截断误差为 ()()()556414.45n nT h y x O h +=+ ()5.12 米尔尼方法也可以通过方程()1.1两端积分()()()()44,n n x n n x y x y x fx y x d x ++-=⎰得到．若将方程()1.1从n x 到2n x +积分，可得()()()()22,.n n x n n x y x y x fx y x d x ++-=⎰ 右端积分通过辛普森求积公式就有()2124.3n n n n n hy y f f f +++=+++ ()5.13 称为辛普森方法．它是隐式二步四阶方法，其局部截断误差为()()()5562.90n nh T y x O h +=-+ ()5.14 7.6方程组和高阶方程7.6.1 一阶方程组前面我们研究了单个方程'y f =的数值解法，只要把y 和f 理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形．考察一阶方程组()'12,,,,i i N y f x y y y = ()1,2,,i N =的初值问题，初始条件给为()00i i y x y = ()1,2,,i N =若采用向量的记号，记()12,,,,TN y y y y =()00012,,,,TN y y y y =()12,,,.TN f f f f =则上述方程组的初值问题可表示为()()'00,,.y f x y y x y ⎫=⎪⎬=⎪⎭()6.1求解这一初值问题的四阶龙格-库塔公式为()1123422,6n n hy y k k k k +=++++ 式中()1,,n n k f x y =21,,22n n h h k f x y k ⎛⎫=++ ⎪⎝⎭32,,22n n h h k f x y k ⎛⎫=++ ⎪⎝⎭()43,.n n k f x h y hk =++或表示为(),11234226i n in i i i i hy y K K K K +=++++ ()1,2,,i N =其中()112,,,,,i i n n n Nn K f x y y y =21112211,,,,,2222i i n n n Nn N h h h h K f x y K y K y K ⎛⎫=++++ ⎪⎝⎭31122222,,,,,2222i i n n n Nn N h h h h K f x y K y K y K ⎛⎫=++++ ⎪⎝⎭41132233,,,,.2i i n n n Nn N h K f x y hK y hK y hK ⎛⎫=++++ ⎪⎝⎭这里in y 是第i 个因变量()i y x 在节点0n x x nh =+的近似值．为了帮助理解这一公式的计算过程，我们再考察两个方程的特殊情形：()()()()''0000,,,,,,.,y f x y z z g x y z y x y z x z ⎧=⎪=⎪⎨=⎪⎪=⎩ 这时四阶龙格-库塔公式具有形式()()112341123422,622.6n n n n h y y K K K K h z z L L L L ++⎫=++++⎪⎪⎬⎪=++++⎪⎭()6.2其中()()()()12113224331211322433,,,,,,222,,,222,,,,,,,,,222,,,222,,.n n n n n n n n n n n nn n n n n n n n n n n n K f x y z h h h K f x y K z L h h h K f x y K z L K f x h y h K z hL L g x y z h hh L g x y K z L h h h L g x y K z L L g x h y h K z h L =⎛⎫=+++ ⎪⎝⎭⎛⎫=+++ ⎪⎝⎭=+++=⎛⎫=+++ ⎪⎝⎭⎛⎫=+++ ⎪⎝⎭=+++⎫⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎭()6.3这是一步法，利用节点n x 上的值,n n y z ，由()6.3式顺序计算1122334,,,,,,,K L K L K L K L，然后代入()6.2式即可求得节点1n x +上的11,n n y z ++． 7.6.2 化高阶方程为一阶方程组关于高阶微分方程（或方程组）的初值问题，原则上总可以归结为一阶方程组来求解．例如，考察下列阶微分方程()()'1,,,,,mm y f x y y y -= ()6.4初始条件为()()()()()11''000000,,,.m m y x y y x y y x y --=== ()6.5只要引进新的变量()1'12,,,,m m y y y y y y -===即可将m 阶方程()6.4化为如下的一阶方程组：()'12'23'1'12,,,,,,,.m m m m y y y y y y y f x y y y -⎫=⎪=⎪⎪⎬⎪=⎪⎪=⎭()6.6初始条件()6.5则相应地化为()()()()1'10020000,,,.m m y x y y x y y x y -=== ()6.7 不难证明初值问题()()6.4,6.5和()()6.6,6.7是彼此等价的．特别的，对于下列二阶方程的初值问题：()()()'''00''00,,,,.y f x y y y x y y x y ⎧=⎪⎪=⎨⎪=⎪⎩引进新的变量'z y =，即可化为下列方程组的初值问题：()()'''00,,,,.y z z f x y z z x y ⎧=⎪=⎨⎪=⎩针对这个问题应用四阶龙格—库塔公式()6.2，有()()112341123422,622.6n n n n h y y K K K K h z z L L L L ++⎧=++++⎪⎪⎨⎪=++++⎪⎩由()6.3式可得1,n K z = ()1,,;n n nL f x y z = 21211,,,;2222n n n n h h h h K z L L f x y K z L ⎛⎫=+=+++ ⎪⎝⎭32322,,,;2222n n n n h h h h K z L L f x y K z L ⎛⎫=+=+++ ⎪⎝⎭()43433,,,.n n n n K z hL L f x h y hK z hL =+=+++如果消去1234,,,K K K K ，则上述格式可表示为()()2112311234,622.6n n n n n h y y hz L L L h z z L L L L ++⎧=++++⎪⎪⎨⎪=++++⎪⎩这里()1,,,n n n L f x y z =21,,,222n n n n h h h L f x y z z L ⎛⎫=+++ ⎪⎝⎭ 2312,,,2242n n n n h h h h L f x y z L z L ⎛⎫=++++ ⎪⎝⎭ 2423,,.2n n n n h L f x h y hz L z hL ⎛⎫=++++ ⎪⎝⎭评注本章研究求解常微分方程初值问题的数值方法．构造数值方法主要有两条途径：基于数值积分的构造方法和基于泰勒展开的构造方法．后一种方法更灵活，也更具有一般性．泰勒展开方法还有一个优点，它在构造差分公式的同时可以得到关于截断误差的估计．基于泰勒展开构造出的四阶龙格-库塔方法（见3.3节）则是计算机上的常用算法．四阶龙格-库塔方法的优点是精度高，程序简单，计算过程稳定，并且易于调节步长．四阶龙格-库塔方法也有不足之处：它要求函数(),f x y 具有较高的光滑性．如果(),f x y 的光滑性差，那么，它的精度可能还不如欧拉公式（2.1节）或改进的欧拉公式（2.4节）．四阶龙格-库塔方法的另一个缺点是计算量比较大，需要耗费较多的机器时间（每一步需四次计算函数(),f x y 的值）．相比之下，汉明方法（5.4节）可以节省计算量（每一步只需两次计算函数(),f x y 的值）．但汉明方法是一种四步法，它不是自开始的，需要借助于四阶龙格-库塔方法提供开始值．对数值方法的分析还涉及到局部截断误差，整体误差，收敛性，相容性及稳定性等概念，特别是有关绝对稳定性的讨论，它涉及计算时步长的选取，如步长选的不合适，舍入误差恶性增长，结果完全错误．本章只对单步法的收敛性和稳。

第七章框架-剪力墙结构在水平荷载下的近似计算方法

第七章框架－剪力墙结构在水平荷载下的近似计算方法本章导学框架：剪力墙结构是由框架和剪力墙组成的一种复合结构体系，它兼具框架结构和剪力墙结构的优点，因而成为高层建筑的主要结构体系。

在水平荷载作用下，因为框架与剪力墙的变形性质不同，不能直接把总水平剪力按抗侧刚度的比例分配到每榀结构上，而是必须采用协同工作方法求得侧移和各自的水平层剪力及内力。

框架剪力墙结构计算的近似方法是将结构分解成平面结构单元，它适用于比较规则的结构，而且只能计算平移时的剪力分配，如果有扭转，要单独进行扭转计算，再将两部分内力叠加。

这种方法概念清楚，结果的规律性较好。

本章主要学习框架：剪力墙结构计算的近似方法，学习中要求同学们熟练掌握协同工作方法的两种计算简图，熟练掌握铰接体系和刚接体系的计算方法的区别与联系。

知识学习第一节概述一．基本假定框剪结构体系在水平荷载作用下的内力分析是一个三维空间超静定问题，通常把它简化为平面结构来计算，并在结构分析中作如下基本假定：①楼板在自身平面内刚度无限大。

这一假定保证楼板将整个计算区段内的框架和剪力墙连成一个整体，在水平荷载作用下，框架和剪力墙之间不产生相对位移。

②当结构体型规则、剪力墙布置比较对称均匀时，结构在水平荷载作用下不计扭转的影响；否则应考虑扭转的影响。

③不考虑剪力墙和框架柱的轴向变形及基础转动的影响。

④结构为线弹性结构。

二．计算简图用连续化解法求总剪力墙与总框架之间的相互作用力，都要解决如何合并总剪力墙、总框架，以及确定总剪力墙和总框架之间的连接和相互作用关系，以便于确定计算简图。

框剪结构用连续化方法求解时，根据连杆刚度情况可以确定两种计算简图：铰接体系和刚接体系。

1.铰接体系在基本假定的前提下，计算区段内结构在水平荷载作用下，处于同一楼面标高处各片剪力墙及框架的水平位移相同。

此时可把平行于水平荷载作用方向的所有剪力墙综合在一起成总剪力墙（一般简化为整体墙），把平行于水平荷载作用方向的所有框架综合在一起成总框架。

数值计算方法第07章数值微分与数值积分

h
2
f '( x) f ( x) f ( x h) f ''( x 2h) h O(h)
h
2
f '( x) f ( x h) f ( x h) 2h
f (3)( x 3h) f (3)( x 3h) h2 O(h2 )
12
心差商公式
sin x2 , cos x2 , sin x , 1 , 1 x3 , ex2 x ln x
17
2. 有些被积函数其原函数虽然可以用初等函数表示,但表达式相当复杂,计算极不方便.
x x1 x0 x1
f
( x0 )
x x0 x1 x0
f
(
x1
)@
x
h
x1
f
( x0 )
x
x0 h
f ( x1 )
则
L1( x)
1 [ h
f
( x0 )
f
( x1 )]
(7.1)
L1( x0 )
1 [ h
f
( x0 )
f
( x1 )],
L1( x1 )
1 [ h
f
( x0 )
f
x0 )( x x2 ) x0 )( x1 x2 )
f
( x1)
(x (x2
x0 )( x x1 ) x0 )( x2 x1 )
f
(x2 )
(x
x1 )( x 2h2
x2 )
f
( x0 )
(x
x0 )( x h2
x2 )
f
(x ( x1 )
x0 )( x 2h2
x1 )
f (x2 )

7 第七章油品静态计量数量计算11.07(最后)

第七章油品静态计量数量计算对油品计量的最终目的是获得其数量（体积或质量）。

物质的质量是由其体积和密度决定的，在油品计量中，油品的密度会由人工或自动采集直接得到，而油品的体积则要通过测得的油品高度查找储油容器容积表得到。

第一节术语及基本计算方法一、术语1、标准温度（t 20）确定某些随温度变化的物理量时选定的一个参照温度，我国规定101.325kPa 大气压下的标准温度为20℃。

2、游离水（FW ）在油品中独立分层并主要存在于油品下面的水，其体积为FW V 。

3、沉淀物和水（SW ）油品中的悬浮沉淀物、溶解水和悬浮水总称为沉淀物和水，其体积为SW V 。

4、总计量体积（to V ）在计量温度下，所有油品、沉淀物和水以及游离水的总测量体积。

5、毛计量体积（go V ）在计量温度下，已扣除游离水的所有油品以及沉淀物和水的总测量体积。

6、毛标准体积（gs V ）在标准温度下，已扣除游离水的所有油品及沉淀物和水的总体积。

7、净标准体积（ns V ）在标准温度下，已扣除游离水及沉淀物和水的所有油品的总体积。

8、表观质量（m ）是油品在空气中称重所获得的数值，也习惯称为商业质量或重量。

它有别于未进行空气浮力影响修正的真空中的质量。

9、毛表观质量（）与毛标准体积对应的表观质量。

10、净表观质量（n m ）与净标准体积对应的表观质量。

11、沉淀物和水的修正系数（CSW ）为扣除油品中的沉淀物和水，将毛标准体积修正到净标准体积或将毛质量修正到净质量的修正系数。

12、罐壁温度修正系数（CTSh ）将油罐从标准温度下的标定容积修正到使用温度下实际容积的修正系数。

13、表观质量换算系数（WCF ）将油品从标准体积换算为空气中的表观质量的系数。

该系数等于标准密度减去空气浮力修正值。

根据国际标准空气浮力修正值为 1.1kg/m 3或0.0011g/cm 3。

即：WCF =20ρ-1.1或WCF =20ρ-0.0011。

计算方法第七章常微分方程数值解法

12
7.1.1 欧拉法及其截断误差
4、欧拉公式的截断误差是O(h2)，公式是1 阶的。
因为
yi+ 1 ? yi
1
h f ( x i , y i ) = y ( x i ) + h y ¢( x i )
1 2
1
2n ) ( n y ( x ) y (( x)i 1 ) ( y i()x i ) y y ( ) ( x h i ) ( )yh ( ( x x i ) xi y x ) ( xi ) x y y 2 n! 2
y0 y( x0 )
i1 i i i
(欧拉公式)
9
7.1.1 欧拉法及其截断误差
例取步长 h=0.1，用欧拉法求解初值问题
ì y ¢= x + y ï ï í ï y (0) = 1 ï î
y i 1 y i h f ( x i , y i ) , i 0 ,1 , 2 , y0 y( x0 )
y f ( x , y ), y( x0 ) y0
x [a , b ]
23
7.1 欧拉法和改进的欧拉法
欧拉公式
y i 1 y i h f ( x i , y i ) , i 0 ,1 , 2 , y0 y( x0 )
h y i 1 y i [ f ( x i , y i ) f ( x i 1 , y i 1 )] , i 0 ,1 , 2 , 2 y y( x ) 0 0
( p)
1.2 ? 1.24
1.528
y 2 = y 1 + 0 .1[( x1 + y 1 ) + ( x 2 + y 2 )] = 1 .2 4 + 0 .1(0 .2 + 1 .2 4 + 0 .4 + 1 .5 2 8) = 1 .5 7 6 8

计算方法引论-第七章

计算方法引论：数值代数⏹解线性方程组的直接法⏹解线性方程组最小二乘问题⏹解线性方程组的迭代法⏹矩阵特征值和特征向量的计算⏹非线性方程及非线性方程组解法第七章线性方程组最小二乘问题•线性最小二乘问题•满秩分解•广义逆矩阵•Gram-Schmidt方法•Householder变换•Givens变换•奇异值分解线性最小二乘问题•线性代数方程组Ax=b(1)–相容:有解, 可能有无穷多解(欠定).–不相容(矛盾,超定):无解.–广义解:最小二乘解.总存在,可能有无穷多.•最小二乘解–求剩余平方和║Ax -b ║2的最小值点–求正规方程(法方程)A T Ax =A T b (通常意义)的解–二者等价:象数据拟合法那样用微分法可得⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++mn mn m m n n n n b x a x a x a b x a x a x a b x a x a x a 22112222212111212111满秩分解与广义逆矩阵•满秩矩阵:A,rank(A)=min{m,n}m×n–行满秩:A,rank(A)=mm×n–列满秩:A,rank(A)=nm×n•满秩分解A= B m×r A r×n,rank(A)=rm×n–(不惟一)可取A的线性无关列为B,它们表出A各列的系数对应为C•广义逆矩阵(惟一)–A+=C T(CC T)-1(B T B)-1B T•注:广义逆矩阵可多个方式定义并确认其惟一性.似乎用奇异值分解更简明实用A+计算•满秩矩阵–行满秩: A+=A T(AA T)-1–列满秩: A+=(A T A)-1A T •非零向量–行向量:x=(x1,x2 ,…,x n)x +=x T/(x12+…+x n2)–列向量x=(x1, x2,…,x n) Tx +=x T/ (x12+…+x n2) •例[][]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=-⨯⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=-⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡--=21211012151112111212211AA +计算又例•A 作满秩分解消元所有列都用1、3列表出131042611713013⎛⎫⎪= ⎪ ⎪⎝⎭A 0.0535710.0178570.0714290.160710.0535710.214290.369050.0119050.380950.422620.0297620.452380.208330.0416670.16667+-⎛⎫⎪- ⎪⎪=-- ⎪- ⎪ ⎪-⎝⎭A 消元矩阵解释•求A +–·–·1–·1131040011100000⎛⎫⎪--⎪ ⎪⎝⎭131040011100111⎛⎫⎪-- ⎪⎪--⎝⎭ ③–②·1 ②–①·2 ③–①·1 1113013210011110⎛⎫⎛⎫⎪=⎪⎪-⎝⎭⎪⎝⎭A 100131042100011111100*********4210011111⎛⎫⎛⎫ ⎪⎪=-- ⎪⎪ ⎪⎪⎝⎭⎝⎭⎛⎫⎛⎫ ⎪= ⎪ ⎪--⎝⎭ ⎪⎝⎭AA+性质•X=A+满足Penrose方程–AXA=A(P1)–XAX=X(P2)–(AX)T=(AX)(P3)–(XA)T=(XA)(P4)•性质–A可逆A+=A-1–(A T)+= (A+)T–(A T A)+=A+(A T)+–(A+A)2=A+A,(AA+)2=AA+•注:不具有逆的某些性质[][]2222))(())((乃知)(1141))(()()()(故2/1)()(,1)(1121)()()(1)(,1121)(,01,11++++++++++++++≠=⎥⎦⎤⎢⎣⎡=≠==⎥⎦⎤⎢⎣⎡==⎥⎦⎤⎢⎣⎡===⎥⎦⎤⎢⎣⎡=xyxyxyxyxyyxyxyxyxxyxyyxyx方程组(1)的解•方程组(1)有解iff AA+b=b–充分性:AA+b=b,则x=A+b满足(1)–必要性:有Ax=b即有AA+(Ax)=b, AA+b=b •(1)有解则其通解为(2)x=A+b+(I-A+A)z, z任意n维向量–(1)有解则A+b是解而A(I-A+A)=A-AA+A=O–设y是解.令z=y-A+b则Az=o.于是z=(I-A+A)z,从而y=A+b+z=A+b+(I-A+A)z.方程组(1)的解(续)•(1)有解时A+b为其通解(2)中惟一2-范数最小者.一般情况下(1)的最小二乘解通解亦(2), A+b仍为其通解中惟一2-范数最小者–(1)有解通解是(2).由于(A+)T(I-A+A)=(A+)T (A+A)T(I-A+A)=(A+)T (A+A-A+A)=O.得║x ║2 = ║A+b║2+║(I-A+A)z║2 +(A+b)T(I-A+A)z= ║A+b║2+║(I-A+A)z║2>0,当(I-A+A)z≠o–一般情况下.令b=c+d, c=AA+b,d=(I-AA+)b,则c T d=0,A T d=o,Ax=c有解y=A+c+(I-A+A)z=A+b+(I-A+A)z,且║b-Ax║2=║c+d-Ax║2 =║c-Ax║2 +║d║2 > ║d║2 ,当Ax≠c.乃证得(1)的最小二乘解通解亦(2).其中A+b为惟一2-范数最小者前己证得.Gram-Schmidt 正交化•G-S 方法可将线性无关的向量组正交化–β1=α1, r 11=║β1║,q 1=β1/r 11–β2=α2-r 12q 1, r 12=(α2, q 1), r 22=║β2║,q 2=β2/r 22–βk =αk -r 1k q 1 -r 2k q 2 -…-r k-1,,k q k , r ik =(αk , q i ), i =1,2, …,k -1, r kk =║βk ║, q k =βk /r kk , k =3, …,n•矩阵表示–A=QR–(α1α2 …αn )=(q 1q 2…q n )⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡nn n n r r r r r r22211111•G-S 方法•修改的G-S 方法–后一算法改变原算法:算出后立即修改使之皆与正交,以后对它们逐个作类似处理.–二算法主要运算量是乘法和加法运算各mn 2次1112/||||=βααfor j = 2:nT T T 112211j j j j j j j --=----βαβαββαββαβ 2/||||j j j =βββ endfor j = 1:n2/||||j j j =βαα for k = j +1:nT k k j k j =-ααβαβ endend1β211,,,n j j j =-：ααααβαβ211,,,2,3,,.n j j j j n =-=：ααααβαβT 211,,,2n j j j j =-=：ααααβαβ1β•A =QR•各列正交化过程5251103202230012---⎛⎫⎪-⎪= ⎪- ⎪-⎝⎭A 0.980580.0377430.176600.0764720.196120.188710.883020.3823600.981310.176600.07647200.397360.91766-⎛⎫⎪-- ⎪= ⎪⎪-⎝⎭Q 5.0990 1.9612 5.49130.588350 2.0381 1.5852 2.528800 2.5166 3.267200.76472---⎛⎫⎪- ⎪=⎪- ⎪⎝⎭R T T 22121(0.076923,0.38462,2,0)=-=-βαβαβ22|||| 2.0381=βT 2(0.037743,0.18871,0.98131,0)=-βT T13235.4913, 1.5852=-=βαβαT T T33131232(0.44444,2.2222,0.44444,1)=--=-βαβαββαβ32|||| 2.5166=βT3(0.17660,0.88302,0.17660,0.39736)=-β-3.2672T4(0.058480,0.29240,0.058480,0.70175)=---β42||||0.76472=βT4(0.076472,0.38236,0.076472,0.91766)=--β12|||| 5.0990=αT1(0.98058,0.19612,0,0)=βT 12 1.9612=-βαT T T 1424340.58835, 2.5288,0.32672=-=-=-βαβαβαHouseholder 变换•定义–H =I -2ww T , ║ w ║ =1•性质–H T =H–H T H =H 2=I–任一x ,║Hx ║=║x ║–任给x 及y ,║y ║=║x ║≠0,总有H 使Hx=y , 不难验证:取w =(x-y )/ ║x-y ║即可.–y 常取坐标轴方向,如y = -sign(x 1 )║x ║ e 1v =x +sign(x 1)║x ║ e 1(w =v / ║v ║)H =I -βvv T ,β=2/v T v用此变换可将矩阵化成上三角(消元)Household变换:算法•变换Hx= -αe1:计算v(存入x)及α=sign(x1)║x║,β–η=max{|xi|}计算v及β时引入的比例因子–xi=x i/η, 1≤i≤n–α= sign(x1)║x║–x1=x1+α, β =(αx1)-1, α= ηα•计算A=HA(H由β,v给出)–设A=(a1… aq)则HA=A-βvv T A=(…a j -βvv T a j…)–算法:对j=1,2,…,qσ=v T aja j =a j -σβv•由此不难导出化上三角的算法Household 正交化•Householder 变换可实现QR 分解–A =QR , Q m ×m 正交阵, R m ×n 上三角阵–实现:作Q p …Q 2Q 1A=R , Q k 是H-变换.p =min{m -1,n }, 即得A =QR , Q =Q 1Q 2…Q p -1•典型步(对照右边矩阵表示)–象消元法那样将右下角矩阵第一列对角元下全变成零(己是则免,H =I )–Ĥ=I -βvv T ,β=2/v T v同前,H =diag(I Ĥ)也是H-变换T⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡*=⎥⎦⎤⎢⎣⎡*⎥⎦⎤⎢⎣⎡v o v o I H A o r A H A O R H I β•Householder 正交化(QR 分解)算法1．输入m n ⨯∈A R ，置1,min(1,).k p m n ==- 2．max(||,,1,,).ik a i k k m η==+ 3．若0η=，则0,0k kk d r ==，否则221/,,,sign()(),1,,,,,,1,,ik ik kk kkmkkk kk k kk kk m j ik ij k i k ij ij j ik a a i k ma a a a a d a r a a d j k n a a a i k m j k n ηαααηαττ-====++=+==-⎛⎫==+ ⎪⎝⎭=-==+∑ 4．若k <p ，则k = k +1，转步骤2；否则，结束. 11213111(1)(2)2223222(1)(2)(3)333333(1)(2)(3)444,,r r r r r r βββ⎛⎫ ⎪⎛⎫⎛⎫ ⎪ ⎪ ⎪=== ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎪⎝⎭A r d v v v v v v v v v (1)11213(1)(2)2223(1)(2)(3)333(1)(2)(3)444,r r r ⎛⎫ ⎪⎛ ⎪==⎪⎪⎝ ⎪⎝⎭A r v v v v v v v v v 运算量:乘、加各次求Q=H 1…H P I 另需乘、加次数各存储方式: 2313mn n -22312()3m n mn n -+2.0198 1.9612 5.49130.588350.20.39223 1.9652 2.157302230012⎛⎫ ⎪- ⎪ ⎪- ⎪-⎝⎭•A =QR•正交化过程k = 1525110320223012---⎛⎫⎪-⎪= ⎪- ⎪-⎝⎭A 5.0990 1.9612 5.49130.588350 2.0381 1.5852 2.528800 2.51663.2672000.76472-⎛⎫⎪-- ⎪= ⎪-⎪-⎝⎭R 0.980580.0377430.176600.0764720.196120.188710.883020.3823600.981310.176600.076472000.397360.91766---⎛⎫⎪-- ⎪= ⎪---⎪-⎝⎭Q a 1 = (5,1,0,0)T ,5η= 1.0198α=T11(2.0198,0.2,0,0)==a v 11110.48548, 5.0990d r β===-a 1 = (1,0.2,0,0)T ,a 11 =1+1.098= 2.019 8用以变换A 的后三列得到1Household 正交化算例(续)•正交化过程k = 2•正交化过程k = 3T2(0.39223,2,0),2η==a 220.19612 1.0190 1.2152a =+=T22(1.2152,1,0)==a v 22220.80755, 2.0381d r β===-2.0198 1.9612 5.49130.588350.2 1.2152 1.5852 2.528801 2.3094 2.6943012⎛⎫⎪- ⎪⎪- ⎪-⎝⎭用以变换A 的最后一列得到用以变换A 的后二列得到T3(2.3094,1), 2.3094η==a T3(1,0.43301)=a 1.0897α=T33(2.0897,0.43301)==a v 33330.43913, 2.5166d r β===-2.01981.9612 5.49130.588350.21.2152 1.58522.528801 2.08973.2672000.433010.76472⎛⎫ ⎪- ⎪ ⎪⎪-⎝⎭T2(0.19612,1,0), 1.0190α==a T 11T 22T330.48548,(2.0198,0.2,0,0)0.80755,(1.2152,1,0)0.43913,(2.0897,0.43301)βββ======v v v R 如前，Q 可由下面的信息生成Givens 变换•定义–G =G (i ,k ,θ)=I +s (e i e k T -e k e i T )+c (e i e i T +e k e k T )•性质–G T G =I–任给x 可使y =G x 的k 分量为零:r =(x i 2+x k 2)1/2 ≠0c =x i /r ,s =x k /r•可用以化上三角形一如消元过程⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-=11),,(c s s c k i G θendend ,11, else ,11, if else 0,1 0 if 22cts tc x x t stc ts x x t x x s c x i k k i ik k =+===+==≥===Givens 正交化算例•算例•过程525110320223012---⎛⎫⎪-⎪= ⎪- ⎪-⎝⎭A 0.980580.0377430.176600.0764720.196120.188710.883020.3823600.981310.176600.076472000.397360.91766--⎛⎫⎪- ⎪=⎪-⎪⎝⎭Q 5.09901.9612 5.49130.5883502.0381 1.5852 2.528800 2.51663.267200.76472---⎛⎫⎪- ⎪=⎪- ⎪-⎝⎭R (2,1)元变为零，t = 0.2, c = 0.980 58, s = 0.196 12.一行：5.099 0 –1.961 2 –5.491 3 –0.588 35二行：0 0.392 23 –1.961 2 2.157 3(3,2)元变零，t = 0.196 12, c = 0.192 45, s = 0.981 31.二行：0 2.038 1 1.585 2 –2.528 8三行：0 0 2.309 4 –2.694 3(4,3)元变零，t = 0.433 01, c = 0.917 66, s = 0.397 36.QR 分解定理•定理设A 是m ×n (m ≥ n ) 矩阵,则A 有QR 分解, 其中Q 是m ×n 的正交矩阵,R 是具有非负对角元的上三角矩阵;而且当m = n 且A 非奇异时R 的对角元皆正上述分解还是唯一的•证⎛⎫= ⎪⎝⎭R A Q O 于是，有T 12T11||||⎛⎫= ⎪ ⎪⎝⎭Q A A v α 1 n –1对(1)(1)m n -⨯-矩阵A 1应用数学归纳法假定，得 212⎛⎫= ⎪⎝⎭R A Q O 其中，Q 2是(m –1)×(m –1)正交矩阵，R 2是具有非负对角元的(1)(1)n n -⨯-上三角矩阵. 这样，令T 121221||||,⎛⎫⎛⎫== ⎪ ⎪ ⎪⎝⎭⎝⎭Q Q R Q R αv 000 则Q 和R 满足定理的要求. 存在性得证.再证唯一性. 设m = n 且A 非奇异，易知R 对角元皆正，假定==A QR QR ，其中，Q , Q 是m m ⨯正交矩阵，R , R 是具有正对角元的上三角矩阵. A 非奇异蕴含着R , R 的对角元均为正数，因此，有 T 1-=Q Q RR 既是正交矩阵又是对角元均为正数的上三角矩阵，只能是单位矩阵. 从而，必有=Q Q ，=R R 即分解是唯一的先证存在性，用数学归纳法. 当n = 1时，定理显然成立. 现假定已经对所有p ×(n –1)矩阵成立，这里假设(1)p n ≥-，设m n ⨯矩阵A 的第一列为1α(可为零向量)，则由定理7.5知，存在m m ⨯正交矩阵Q 1，使得T 11121||||=Q e αα于是，有T 12T11||||⎛⎫=⎪ ⎪⎝⎭0Q A A v α 1 n –1 对(1)(1)m n -⨯-矩阵A 1应用数学归纳法假定，得212⎛⎫= ⎪⎝⎭R A Q O 1m –1最小二乘解:列满秩•列满秩时求(1)的最小二乘解–形成正规方程A T Ax=A T b(n阶)(乘法和加法各mn2/2次) 用平方根法(乘法和加法各n3/6次)用G-S作A=QR:R T Rx=R T Q T b,Rx=Q T b(各mn2次) –用Householder变换或Givens变换作QR分解║Ax-b║2 =║Q T Ax-Q T b║2==║Rx-c1║2+║c2║2,解x=R-1c1,最小值║c2║2注: 若记Qm×m =(Q1 Q2), Q1是m×n阵则有c1=Q1b,c2=Q2b⎥⎦⎤⎢⎣⎡-⎥⎦⎤⎢⎣⎡21ccxOR最小二乘解:列满秩算例•求最小二乘解•求最小二乘解L123525110320223012x x x ---⎛⎫⎛⎫⎛⎫ ⎪ ⎪- ⎪ ⎪ ⎪= ⎪ ⎪ ⎪- ⎪ ⎪ ⎪⎝⎭-⎝⎭⎝⎭12326102831081442814399x x x ---⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪-=- ⎪ ⎪ ⎪ ⎪ ⎪ ⎪--⎝⎭⎝⎭⎝⎭正规方程5.0190001.9912 2.038105.4913 1.58522.5166⎛⎫⎪- ⎪⎪-⎝⎭T( 1.6023,0.23099, 1.2982)=---x 5.09901.9612 5.49130.5883502.0381 1.5852 2.528800 2.51663.26720000.76472-⎛⎫ ⎪-- ⎪=⎪- ⎪-⎝⎭R 5251103202230012---⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭A 方程H-正交化增广矩阵T ( 1.6023,0.23099, 1.2982)=---x 剩余是0.764 721.9612奇异值分解(SVD)•矩阵奇异值–A m ×n 的奇异值σ1≥σ2≥…≥σr >σr+1=…=σn =0是A T A 的特征值λ1, …,λn 的平方根•奇异值分解定理),,diag(正交阵,1Tr n n m m σσΣV U V O O O ΣU A =⎥⎦⎤⎢⎣⎡=⨯⨯形式阵即得矩,,1,连同,得,,,1由),,,(阵正交增.则,,1,/令,,),,,(证.取TTT21T TT21分解補成r k u Av o Av o Av A v o Av A n r k u u u U u u r k Av u v Av A I V V v v v V k k k k k kk m kj j k k k k k k k n =====+========σδσλ奇异值分解(续)•推论(记号同前)–分解形式A=σ1u1v1T+σ2u2v2T+…+σr u r v r TA=ÛΣŴT, Û=(u1u2…u r),Ŵ =(v1v2…v r)–空间关系R(A)=Span{u1,u2,…,u r}N(A)=Span{v r+1,v r+2,…,v n}R(A T)=Span{v1,v2,…,v r}N(A T)=Span{u r+1,u r+2,…,u m}R(A)=N(A T)⊥, R(A T)=N(A)⊥SVD 与A+•X =A +满足Penrose 方程–AXA =A (P1)–XAX =X (P2)–(AX )T =(AX )(P3)–(XA )T =(XA )(P4)•由SVD 解出X =A +T T 111T111rr ru v u v UO O O ΣV A σσ++=⎥⎦⎤⎢⎣⎡=-+T 11T:乃得唯一解得(P4)代入得(P3)代入得(P2)代入得(P1)代入对应分块SVD 解.取U O OO ΣV X O L O K K L ΣM ΣS U M L K S V X ⎥⎦⎤⎢⎣⎡=====⎥⎦⎤⎢⎣⎡=--图示A与A+•SVD绐出A与A+在标准正交基下向量对应关系V1= Span{v,v2,…,v r}, V2= Span{v r+1,v r+2,…,v n},1U1= Span{u,u2,…,u r}, U2=Span{u r+1,u r+2,…,u m}1最小二乘解•(1)的最小二乘解–推导:按前述标准正交基分解,再求║Ax-b║2最小x=c1v1+c2v2+…+c n v nA x=c1σ1u1+c2σ2u2+…+c rσr u rb=u1T b u1+u2T b u2+…+u m T bu m║Ax-b║2最小:c=u k T b/σk,k≤r,其余任意k–最小二乘解通解x=1/σ1u1T b v1+1/σ2u2T b v2+…+1/σr u r T b v r+ v r +…+c n v n, c r+1, …, c n任意c r+1–最小2-范数最小二乘解y=1/σ1u1T b v1+1/σ2u2T b v2+…+1/σr u r T b v rSVD与最小二乘解•上述结果亦可借助SVD得到–A, A+代入通解(2)x=A+b+(I-A+A)z(z任意n维向量)2-范数最小A+b=(1/σ1v1u1T+…+1/σr v r u r T)b(I-A+A)z=(I-(v1v1T+v2v2T+…+v r v r T) )z=(v r+1v r+1T+…+v n v n T)z =c r+1v r+…+c n v n –由SVD直接推出最小二乘解║Ax-b║2=║UΛV T x-b║2=║ΛV T x-U T b║2=║Λc-U T b║2 ,这里Λ=diag(ΣO),c=V T x的i,U T b的i分量u i T b.从而可得结果.分量ci代入正规方程A T Ax=A T b关于A +的定义•A +有多个等价定义–由满秩分解:C T (CC T )-1(B T B )-1B T–由Penrose 方程.–由SVD:–由最小二乘解:(1)中任一b 对应唯一最小2-范数最小二乘解x 所确定的矩阵.–由线性算子确定的矩阵.线性算子f :R m →R nf (y )=x ,当y ∈R (A ), Ax =y f (y )=o ,当y ⊥R (A )•注.一个或几个Penrose 方程可定义多种广义逆T1V O O O ΣU ⎥⎦⎤⎢⎣⎡-。

微积分的数值计算方法数值微分

r=dx/dt/x =0.0283 0.0146 0.0152 0.0097 0.0093 0.0166 0.0154 0.0113 0.0100 0.0113
将节点处的增长率作三次样条插值
年份增长率 1900 0.0283 1901 0.0255 1902 0.0230 1935 0.0082 1936 0.0081 1937 0.0083 1953 0.0172 1954 0.0172 1979 0.0100 1980 0.0100 1981 0.0109 1989 0.0111 1990 0.0113
f ( x 0 ) 21h(3f04f1f2) f ( x n ) 21h(fn24fn13fn)
--------(11)
称(11)式为分段三点公式
实际中下面的公式很有用
f
(
xk
xk1 2
)
1( h
f k 1
fk
)
例: 回到实例(美国人口)
1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4
E 2(x0)f(3 3 )!()(x0x1)x (0x2)
ห้องสมุดไป่ตู้h2 3
f (3)( )
E 2(x1)f(3 3 )!()(x1x0)x (1x2)
h2 6
f (3)()
E 2(x2)f(3 3 )!()(x2x0)x (2x1)
h2 3
f (3)( )
f ( x0 )
21h(3f04f1f2)
1( h
f1
f0 )
h f (2)( )

计算方法第7章习题 - 参考答案

参考答案：习题七7.1 运用Euler 方法和改进Euler 方法求下列初值问题在给定区间上的数值解, 计算结果保留四位小数。

(1) ⎪⎩⎪⎨⎧=≤≤=-=04.0,2.00,1)0(22h x y y x dxdy； (2) ⎪⎩⎪⎨⎧=≤≤=-=1.0,5.00,1)0(h x y ydxdy。

解：(1) 5,4,3,2,1,0,,04.0,1)0(,),(22====-=n nh x h y y x y x f n8360.08635.08935.09262.09616.01Euler 8299.08583.08894.09232.096.01Euler 2.016.012.008.004.00改进k x (2) 5,4,3,2,1,0,,1.0,1)0(,),(====-=n nh x h y y y x f n6071.06708.07412.0819.0905.01Euler 5905.06561.0729.081.09.01Euler 5.04.03.02.01.00改进k x 7.2 用Euler 方法和改进Euler 方法求初值问题⎪⎩⎪⎨⎧=+=0)0(y bax dx dy的解在),2,1(, ==n hn x n 处的近似值。

bnh ah n n bnh n ah bh anh h n a y bhanh y b anh h y b ax h y y x hf y y y x b ax y x f n n n n n n n n n ++=+++++==++-+=++=++=++=+===+=-+2222121002)1()...210(...2)1()()(),(0,0,),(7.3 运用标准四阶Runge--Kutta 法求初值问题⎪⎩⎪⎨⎧=+=1)1(32y yx x dx dy的解在x =1.1,1.2,1.3处的近似值, 计算结果保留三位小数。

102.2)3.1(,587.1)2.1(,24.1)1.1(===y y y7.4 运用标准四阶Runge--Kutta 法求初值问题⎪⎩⎪⎨⎧=--=1)0(2y xy y dx dy在区间[0,1]上的数值解, 取步长h =0.2, 将计算结果与准确值1)12()(---=x e x y x 进行比较。

第7章框架-剪力墙结构近似计算方法

EI
梁刚域的长度，墙宽为 hc
hc hb hc hb hc hb hc hb al = − ， bl = − ，a = ( − )/l ，b = ( − )/l 。 2 4 2 4 2 4 2 4
1. 两端有刚域的连梁刚度 mab =
(1 + a − b) 6 EI (1 − a + b) 6 EI ， mba = 。 3 3 (1 + β )(1 − a − b) l (1 + β )(1 − a − b) l
图 7-2 框架-剪力墙结构铰接体系
当与水平力作用方向相同的剪力墙在自身平面内通过连梁与框架或其他同一方向剪力墙相联结时，剪力墙与连梁的结点为刚接点，除此之外均为铰接点。
图 7-3 框架-剪力墙结构刚接体系
方向水平力作用方向下（图7-4（a）所示结构平面，在y方向水平力作用方向下（图7-4（b）），（）所示结构平面，（）），可以将所有纵向剪力墙和纵梁全部删除后，方向有4片剪力墙片剪力墙。个框架中可以将所有纵向剪力墙和纵梁全部删除后， y方向有片剪力墙。8个框架中其在两侧有框架梁。个框架边柱其只在一侧有框架梁。个框架边柱，柱，其在两侧有框架梁。8个框架边柱，其只在一侧有框架梁。总连梁包含 4个刚结点。个刚结点。个刚结点方向水平力作用方向下（）），可以将所有横向剪力墙和横在x方向水平力作用方向下（图7-4（c）），可以将所有横向剪力墙和横方向水平力作用方向下（）），梁全部删除后，片剪力墙。个框架中柱其在两侧有框架梁。个框个框架中柱，梁全部删除后，有2片剪力墙。14个框架中柱，其在两侧有框架梁。6个框片剪力墙架边柱，其只在一侧有框架梁。总连梁包含2个刚结点个刚结点。架边柱，其只在一侧有框架梁。总连梁包含个刚结点。

三重积分概念及部分计算方法

D
D z1 ( x, y)
4
z
若在 xoy 平面上的投影区域记为Dxy，则有
f ( x, y, z)dv dxdy z2 ( x, y) f ( x, y, z)dz

Dxy
z1 ( x, y)
o
a
投影区域Dxy用不等式表示:
b x
1(x) y 2(x)，a x b
dy
2 ( x) dy
z2 ( x, y) f ( x, y, z)dz

a
1 ( x)
z1( x,y)
5
f ( x, y, z)
b
dy
2 ( x) dy
z2 ( x, y) f ( x, y, z)dz

a
1 ( x)
z1( x,y)
公式(2)把三重积分化为先对z、次对y、最后对x的三次积分。
y
则将二重积分化为二次积分，于是得到三重积分的计算公式：
S2 : z z2(x, y) z2

z1 S1 : z z1( x, y)
D (x, y)
y y1( x)
y
y y2(x)
y 2(x)
Dxy y 1(x)
oa
bx
f ( x, y, z)
b
－∞< z < +∞
o

②三组坐标面分别为
x
y
r ＝常数，即以z 轴为轴的圆柱面；
M (r, , z)
• M(x,y,z)
z
r
y
x
•P(r, ,0)
＝常数，即过z 轴的半平面；