非线性方程组的解法
合集下载
数值分析lec非线性方程组解法
第27页/共36页
取初始值（1，1，1），计算如
下N
x
y
z
0 1.0000000 1.0000000 1.00000000
1 2.1893260 1.5984751 1.3939006
2 1.8505896 1.4442514 1.2782240
3 1.7801611 1.4244359 1.2392924
)
第17页/共36页
第18页/共36页
简单迭代法
将方程组 F(X)=0, 写成与之等价的形式:X=G(X), 然后再 利用 X(k+1)=G(X(k)), k=1,2,…,求解原方程的根。
例子：
4
x1
x2
0.1ex1
1
x1
4x21 8x120x1(k
1)
1 (1 4
x(k) 2
3x2(k ) 12x2(k
5, ) 5
.
逐次迭代得结果.
第26页/共36页
例：用牛顿法解方程组
xy z2 1 xyz y2 x2 2
ex z ey 3
xy z2 1
F
(x)
xyz
y2
x2
2
ex z ey 3
y
F
( x)
yz
ex
x
2 z
xz 2 y
xz
e y
1
来近似f(x)在xk处的一阶导数
xk1
xk
f (xk )(xk x0 ) f (xk ) f (x0 )
第11页/共36页
牛顿下山法：
目的是解决初值的选取范围太小这一困难。
构 造 迭 代 格 式 为 ：xk 1
第十章非线性方程及非线性方程组解法
(
x
)
n
lim
n
x
n
若
{x
}
n
收敛，即
lim xn x*，则：
n
x* (x*) f (x*) 0
迭代过程的几何表示
x (x) :
y x 交点即真根。
y (x)
yx
y
Q1
Q2
P* P2
O x* x2
P1
x1
y (x)
P0
x0
x
例：求方程 f (x) x3 x 1 0 在x0 1.5附近的根x*. 解：（1） 将方程改写为 x 3 x 1
第十章 方程求根
求解非线性方程
f (x) 0 f 是非线性函数，
例：代数方程
a x a x a x a f (x) n
n1 L
0, n 1。
n
n1
1
0
例: 超越方程
f (x) ex sin x 0
§1. 非线性方程实根的对分法(二分法)
设 f (x) 在[a,b] 上连续且 [a,b] 有且仅有一个根又
xn1 (xn ) (n 0,1,L )
均收敛于x*，并有
x* xn
Ln 1 L
x1 x0
收敛充分性定理(一、2)
证：由条件（2）知(x)在[a, b]上连续。 令 (x) x (x),则 (x)在[a,b]上连续，且
(a) a (a) 0, (b) b (b) 0 故存在 [a,b]，使得（） 0,即 （）, 所以方程x (x)在[a,b]内有根。
可先用二分法或经验确定迭代初值x0 0.5，再按牛
顿公式进行迭代。
Newton法具有收敛快，稳定性好，精度高等优点，是求 解非线性方程的有效方法之一。但它每次迭代均需计算函 数值与导数值，故计算量较大。而且当导数值提供有困难 时， Newton法无法进行。
第7章非线性方程组的数值解法
( 1, 1 )
f 1 y f 2 2 y
2 y ( 1,1 ) 2
( 1,1 )
( y 3) ( 1, 1 )
( 1, 1 )
( x 1) ( 1 , 1 ) 2
( 1,1 )
f 1 f 2 2 2[ 2 * ( 3) ( 2 ) * ( 2 )] 4 f1 f2 g10 x ( 1,1) x ( 1,1) x f 1 f 2 g 2 2[ 2 * ( 3) 2 * ( 2 )] 20 20 y y f 1 y f 2 ( 1, 1 ) ( 1, 1 )
完
f ( x0 h, y0 k ) f ( x0 , y0 ) ( h k ) f ( x0 , y0 ) x y 1 2 ( h k ) f ( x 0 , y0 ) 2! x y 1 n ( h k ) f ( x 0 , y0 ) n! x y 1 n 1 ( h k ) f ( x0 h, y0 k ) ( n 1)! x y
2
2
令
0
得 f 1 f 1 ( g10 x g 20 y ) f 1 ( g10 ( g f 1 g f 1 ) 2 ( g 10 20 10 x y f 2 g 20 x f 2 g 20 x f 2 ) f2 y f 2 2 ) ( x y
1
f 1 ( x 0 , y0 ) f ( x , y ) 2 0 0
从n到n+1的迭代格式为：
f 1 ( x n , y n ) xn 1 x n x y y f 2 ( xn , yn ) n 1 n x
f 1 y f 2 2 y
2 y ( 1,1 ) 2
( 1,1 )
( y 3) ( 1, 1 )
( 1, 1 )
( x 1) ( 1 , 1 ) 2
( 1,1 )
f 1 f 2 2 2[ 2 * ( 3) ( 2 ) * ( 2 )] 4 f1 f2 g10 x ( 1,1) x ( 1,1) x f 1 f 2 g 2 2[ 2 * ( 3) 2 * ( 2 )] 20 20 y y f 1 y f 2 ( 1, 1 ) ( 1, 1 )
完
f ( x0 h, y0 k ) f ( x0 , y0 ) ( h k ) f ( x0 , y0 ) x y 1 2 ( h k ) f ( x 0 , y0 ) 2! x y 1 n ( h k ) f ( x 0 , y0 ) n! x y 1 n 1 ( h k ) f ( x0 h, y0 k ) ( n 1)! x y
2
2
令
0
得 f 1 f 1 ( g10 x g 20 y ) f 1 ( g10 ( g f 1 g f 1 ) 2 ( g 10 20 10 x y f 2 g 20 x f 2 g 20 x f 2 ) f2 y f 2 2 ) ( x y
1
f 1 ( x 0 , y0 ) f ( x , y ) 2 0 0
从n到n+1的迭代格式为：
f 1 ( x n , y n ) xn 1 x n x y y f 2 ( xn , yn ) n 1 n x
非线性代数方程组的解法
δ 1 = δ 0 + Δδ 1
11
(K1 )−1
=
Δδ 1 F1 − F0
= δ1 −δ0 ψ 1 −ψ 0
Δδ 2 = (K1)−1(R − F1)
………
(K i )−1 =
Δδ i Δψ i
=
δ ψ
i i
− δ i−1 −ψ i−1
(2.13)
显然 K i 就是相应于 Δδ i = δ i − δ i−1 与 Δψ i =ψ i −ψ i−1 的割线劲度。但实际上对于多维情
是对称矩阵。所以式（2.21）是对称秩 1 算法。 (2) 秩 2 算法
一个 N × N 阶的秩 2 矩阵，总可以表示为
[ ] 1T B2T
⎤ ⎥ ⎦
=
A1 B1T
+
A2 B2T
(2.22)
式中A1、A2、B1和B2均为N×1 维向量。将上式代入（2.14），再代入（2.15）得 A1B1T Δψi + A2 B2T Δψi = Δδi − (K i−1 )−1 Δψi
ψ(δi ) ≡ K (δi )δi − R ≠ 0
ψ(δ) 作为对平衡偏离的一种度量，称为失衡力。
对于一个单变量问题的非线性方程，直接迭代法的计算过程如图 2.1 和图 2.2 所示，它
们分别给出 F~δ为凸和凹曲线时的迭代过程。可以看出 K(δ)就是过曲线上点(δ, F(δ)与原点
的割线斜率。对于单变量问题，这一迭代过程是收敛的，但对多自由度情况，由于未知量通
似解。若
ψi−1 = ψ(δi−1) ≡ F (δi−1) − R ≠ 0
希望能找到一个更好的、方程（2.4）的近似解为
δ = δi = δi−1 + Δδi
非线性方程的解法
20世纪60年代中期以后，发展了两种求解非线性方程组(1)的新方法。
一种称为区间迭代法或称区间牛顿法，它用区间变量代替点变量进行区间迭代，每迭代一步都可判断在所给区间解的存在惟一性或者是无解。
这是区间迭代法的主要优点，其缺点是计算量大。
另一种方法称为不动点算法或称单纯形法，它对求解域进行单纯形剖分，对剖分的顶点给一种恰当标号，并用一种有规则的搜索方法找到全标号单纯形,从而得到方程(1)的近似解。
这种方法优点是，不要求f(□)的导数存在,也不用求逆，且具有大范围收敛性，缺点是计算量大编辑摘要目录• 1 正文• 2 牛顿法及其变形• 3 割线法• 4 布朗方法• 5 拟牛顿法•非线性方程组数值解法 - 正文n个变量n个方程(n >1)的方程组表示为(1)式中ƒi(x1,x2,…,x n)是定义在n维欧氏空间R n的开域D上的实函数。
若ƒi中至少有一个非线性函数，则称(1)为非线性方程组。
在R n中记ƒ＝则(1)简写为ƒ(尣)=0。
若存在尣*∈D，使ƒ(尣*)＝0，则称尣*为非线性方程组的解。
方程组(1)可能有一个解或多个解，也可能有无穷多解或无解。
对非线性方程组解的存在性的研究远不如线性方程组那样成熟，现有的解法也不象线性方程组那样有效。
除极特殊的方程外，一般不能用直接方法求得精确解，目前主要采用迭代法求近似解。
根据不同思想构造收敛于解尣*的迭代序列{尣k}(k=0，1，…)，即可得到求解非线性方程组的各种迭代法，其中最著名的是牛顿法。
非线性方程组数值解法 - 牛顿法及其变形牛顿法基本思想是将非线性问题逐步线性化而形成如下迭代程序：(2)式中是ƒ(尣k)的雅可比矩阵,尣0是方程(1)的解尣*的初始近似。
这个程序至少具有2阶收敛速度。
由尣k算到尣k+的步骤为:①由尣k算出ƒ(尣k)及;②用直接法求线性方程组的解Δ尣k；③求。
由此看到迭代一次需计算n个分量函数值和n2个分量偏导数值，并求解一次n阶线性方程组。
非线性方程组的数值解法及最优化方法课件
拟牛顿法求解非线性方程组
拟牛顿法是牛顿法的改进，通过构造一个近似于真实Hessian矩阵的对称正定矩阵来逼近， 从而加快了算法的收敛速度。
信赖域方法求解非线性方程组
信赖域方法是一种基于梯度信息的迭代算法，通过在每一步中计算一个小的搜索方向，并 限制步长，以避免算法发散。
最优化方法案例
梯度下降法求解无约束最优化问题
梯度下降法是一种迭代算法，通过不断沿负梯度方向更新变量，最终找到最优化问题的最小值点。该方法适用于求解 无约束最优化问题。
牛顿法求解无约束最优化问题
牛顿法是一种基于二阶导数的迭代算法，通过不断逼近函数的极小值点，最终求解无约束最优化问题。该方法适用于 求解具有多个局部最小值的问题。
遗传算法求解约束最优化问题 遗传算法是一种基于生物进化原理的随机搜索算法，通过模拟生物进化过程中的自然选择和遗传机制， 在解空间中进行高效搜索，最终找到满足约束的最优解。
和稳定性。
约束最优化方法
拉格朗日乘数法
通过引入拉格朗日函数，将约束最优化问题转化为无 约束最优化问题求解。
罚函数法
通过引入罚函数，将约束条件转化为无约束条件，通 过迭代更新求解。
序列二次规划法
结合拉格朗日乘数法和牛顿法的思想，通过迭代逼近 最优解。
混合整数最优化方法
01
02
03
分支定界法
将整数约束转化为区间约 束，通过不断分支和剪枝 来逼近最优解。
非线性方程组与最优化方法的结合案例
非线性规划问题
非线性规划是最优化领域中一类重要的数学问题，其目标函数和约束条件都是非线性的。常见的非线性规划问题 包括最小二乘问题、二次规划问题等。求解非线性规划问题的常用方法包括梯度下降法、牛顿法等。
拟牛顿法是牛顿法的改进，通过构造一个近似于真实Hessian矩阵的对称正定矩阵来逼近， 从而加快了算法的收敛速度。
信赖域方法求解非线性方程组
信赖域方法是一种基于梯度信息的迭代算法，通过在每一步中计算一个小的搜索方向，并 限制步长，以避免算法发散。
最优化方法案例
梯度下降法求解无约束最优化问题
梯度下降法是一种迭代算法，通过不断沿负梯度方向更新变量，最终找到最优化问题的最小值点。该方法适用于求解 无约束最优化问题。
牛顿法求解无约束最优化问题
牛顿法是一种基于二阶导数的迭代算法，通过不断逼近函数的极小值点，最终求解无约束最优化问题。该方法适用于 求解具有多个局部最小值的问题。
遗传算法求解约束最优化问题 遗传算法是一种基于生物进化原理的随机搜索算法，通过模拟生物进化过程中的自然选择和遗传机制， 在解空间中进行高效搜索，最终找到满足约束的最优解。
和稳定性。
约束最优化方法
拉格朗日乘数法
通过引入拉格朗日函数，将约束最优化问题转化为无 约束最优化问题求解。
罚函数法
通过引入罚函数，将约束条件转化为无约束条件，通 过迭代更新求解。
序列二次规划法
结合拉格朗日乘数法和牛顿法的思想，通过迭代逼近 最优解。
混合整数最优化方法
01
02
03
分支定界法
将整数约束转化为区间约 束，通过不断分支和剪枝 来逼近最优解。
非线性方程组与最优化方法的结合案例
非线性规划问题
非线性规划是最优化领域中一类重要的数学问题，其目标函数和约束条件都是非线性的。常见的非线性规划问题 包括最小二乘问题、二次规划问题等。求解非线性规划问题的常用方法包括梯度下降法、牛顿法等。
42 非线性方程组的迭代解法讲解
* (k )
x ( k ) x ( k 1) x
(k )
；
2o 由
L知简单迭代法是线性收敛的；
3o 对线性方程组迭代函数G ( x ) Bx d , 有L= B <1是收敛的充分 必要条件。
局部收敛定理 定理5(局部收敛定理 ) 设G：D R n R n ,x * int( D )
其中, 0 k 1, k 1, 2,
, n。
三、收敛向量序列的收敛速度
定义3 设向量序列 xk 收敛于 x * , ek x * xk 0,
k 1,2,
, 如果存在常数r 1和常数c 0，使极限
lim
k
e
k
e k 1
r
c
r
成立，或者使得当k K (某个常数)时，有 ek 1 ek
（4Байду номын сангаас2.2）
其中，F : D R n R n是定义在区域D R n上的向量 值函数。 若存在x * D , 使F ( x * ) ，则称x *是方程组(4.2.1)或 (4.2.2)的解。
二、多元微分学补充
定义1 设f ：D R n R，x int( D ) (即x是D的内点)， 若存在向量l ( x ) R n ,使极限
L (k ) ( k 1) L(1 L ) ( k ) ( k 1) x x x x 1 L 1 L L * (k ) 再让m , 得 x x x ( k ) x ( k 1) ■ 1 L
m
i 1 i 1
说明
1o 简单迭代法的精度控制与终止条件e( k ) x * x ( k +1) x x
x ( k ) x ( k 1) x
(k )
；
2o 由
L知简单迭代法是线性收敛的；
3o 对线性方程组迭代函数G ( x ) Bx d , 有L= B <1是收敛的充分 必要条件。
局部收敛定理 定理5(局部收敛定理 ) 设G：D R n R n ,x * int( D )
其中, 0 k 1, k 1, 2,
, n。
三、收敛向量序列的收敛速度
定义3 设向量序列 xk 收敛于 x * , ek x * xk 0,
k 1,2,
, 如果存在常数r 1和常数c 0，使极限
lim
k
e
k
e k 1
r
c
r
成立，或者使得当k K (某个常数)时，有 ek 1 ek
（4Байду номын сангаас2.2）
其中，F : D R n R n是定义在区域D R n上的向量 值函数。 若存在x * D , 使F ( x * ) ，则称x *是方程组(4.2.1)或 (4.2.2)的解。
二、多元微分学补充
定义1 设f ：D R n R，x int( D ) (即x是D的内点)， 若存在向量l ( x ) R n ,使极限
L (k ) ( k 1) L(1 L ) ( k ) ( k 1) x x x x 1 L 1 L L * (k ) 再让m , 得 x x x ( k ) x ( k 1) ■ 1 L
m
i 1 i 1
说明
1o 简单迭代法的精度控制与终止条件e( k ) x * x ( k +1) x x
5-非线性方程组的数值解法及最优化方法
1 4 0
非线性方程组的数值解法
x10=0; x20=0; k=0; while 1 k=k+1; x1k=(1+x20-0.1*exp(x10))/4; x2k=(x10-x10^2/8)/4; %雅克比迭代法 %x2k=(x1k-x1k^2/8)/4; %高斯-赛德尔迭代法 err1=abs(x1k-x10); err2=abs(x2k-x20); err=max(err1,err2); if err<=0.00000000005 break; end x10=x1k; x20=x2k; end
0.0000055305 0.0000001511 0.0000000041 0.0000000001
非线性方程组的数值解法
牛顿迭代法：根据求解非线性方程的牛顿迭代法，如果已经 k k T ，则 ,, xn 给出方程组 Fx 0 的一个近似根 xk x1k , x2 可把函数 Fx 的分量 fi x, i 1,2,, n 在 x k 处按多元函数泰 勒公式展开，取其线性部分做近似，得
(0.2325668498,0.0564514831) (0.2325670008,0.0564515487) (0.2325670050,0.0564515196) (0.2325670051,0.0564515197) (0.2325670051,0.0564515197)
0.0002023950
所以有
1 x φx 1 2 x1
0
T
取初值 x 代公式收敛。
T 0 x 0 , 0 附近 φx 1，所以迭 0,0 ，在
1 1 x 1 e 40 x2 2 1 1 x1 x2 4 16
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
非线性方程组的解法
非线性方程组的解法
主要有两类： 几何非线性，材料非线性 几何非线性 物理位置：平衡方程必须按照变形后的几何位置建立 数学描述：物体的变形一般由位移的一阶微分求得， 当变形很大，高阶微分不能忽略 举例：
非线性方程组的解法
主要有两类： 几何非线性，材料非线性 材料非线性 混凝土，钢筋，粘结滑移，时变效应 一般认为，混凝土开裂即存在材料非线性 其他非线性：边界非线性 常见的是接触问题
[K ][∆δ ] = [∆P ]
n n n
[∆δ n ] = [K n ]−1 [∆Pn ]
[δ n ] = [δ n−1 ] + [∆δ n ]
步骤：……. 缺点：存储两个刚度，存贮量大
非线性方程组的解法
逐步增量法－中点刚度增量法
施加一半荷载，求出相应的位移（中点位移）
[K n−1 ]∆δ n − 1 = 1 [∆Pn ]
非线性方程组的解法
逐步搜索法 可求出极限荷载，无法求解
步骤：1. 加一级荷载 ∆P ，结果发散，退回 2. 加 1 ∆P ，若结果发散，再次退回
2
P −δ
下降段
3. 加
1 ∆P ，……至收敛，得到极限荷载 4
非线性方程组的解法
虚加刚性弹簧法 加虚拟弹簧, 改变原结构形式，新的结构形式 曲线在一定范围内没有下降段
[B] = [M ][R][Q]
其他数值方法
刚体弹簧元
k n [D] = 0
Tresca屈服条件
0 ks
本构矩阵 (2x2)
平面应变下的莫尔-库仑屈服条件
平面应力下的莫尔-库仑屈服条件
其他数值方法
刚体弹簧元
其他数值方法
无网格法
特点 在计算域上，用离散的点拟合场函数，采用基于节点（非单元节点） 的位移插值函数 部分或彻底地消除网格，避免了有限单元法中单元网格的限制 在处理大变形、开裂等问题时，有其优势 采用离散的点拟合场函数的方法 核（权）函数法 移动最小二乘法 径向基函数法 然后建立方程求解
[δ 1 ] = [K 0 ]−1 [P]
节点荷载 [P1 ] ， 即时切线刚度 [K 1 ]
[∆P1 ] = [P] − [P1 ]
[∆δ 2 ] = [K1 ]−1 [∆P1 ]
[δ 2 ] = [δ1 ] + [∆δ 2 ]
[δ n ] = [δ n−1 ] + [∆δ n ]
节点荷载 [P2 ] ， 即时切线刚度 [K 2 ]
非线性方程组的解法
添加行列法
非线性方程组的解法
强制迭代法 结构总刚是负定的,只要施加负荷载,可求得相应的位移增量
问题:适用范围?
非线性方程组的解法
硬化刚度法 人为增大结构刚度,改善收敛性 一般结构切线刚度可分为
KT = K E − K p
弹性 非线性影响 刚度 刚度
K T = K E − λK p
能量准则
([P ] [δ ]) res
T
两曲线的接近如何描述？
非线性方程组的解法
两个问题： 考虑结构负刚度的一些算法 弹塑性单元应力调整计算
非线性方程组的解法
考虑结构负刚度的一些算法 问题：混凝土结构 P − δ 下降段，由力求位移求解困难 原因：上升段，刚度矩阵正定，线性方程组有解 下降段，刚度矩阵非正定，求解困难 负刚度：非正定的刚度矩阵 几种方法：逐步搜索法，虚加刚性弹簧法，位移加载法， 强制迭代法，硬化刚度法，弧长法
刚体弹簧元 刚体＋弹簧 自由度少，总刚阶数少 可在在极限分析领域运用 思路 弹簧单元刚度矩阵 界面应力应变关系 界面法向和切向的应力应变关系 界面法向和切向刚度
其他数值方法
刚体弹簧元
[K ] = ∫ [B ] [D ][B ]ds
T
k n [D] = 0
0 ks
单元刚度矩阵 (6x6) 本构矩阵 (2x2) 形函数矩阵 (2x6)
网格大小与数目的选择 块体大小，计算精度
其他数值方法
公共面法 块体间接触判断 公共面的定义（如何得到） 假象面 令两块体向薄板靠近 两块体在接触时停止运动 薄板被固定在一个位置，该面即为公共面 块体与公共面间的接触关系
其他数值方法
公共面法 接触类型
其他数值方法
颗粒离散元法 多边形、多面体 圆形、球形 可以以链杆或弹簧连接颗粒 质点－桁架模型
2
2
δ 1 = [δ n −1 ] + ∆δ 1 n− 2 n− 2
K n− 1 确定相应的刚度（中点刚度） 2
以该刚度作为初始刚度求当前位移增量
K n − 1 [∆δ n ] = [∆Pn ] 2
优点: 比两步欧拉折线法精度高
P −δ
问题:在何处加弹簧? 弹簧刚度如何?
非线性方程组的解法
位移控制法 支座位移控制加载,可求支座的反力, 和全部 P − δ 曲线 两个算法:刚度矩阵重新排列法, 添加行列法
非线性方程组的解法
添加行列法
[p1 p2]T为参考矢量, ∆λ 为荷载控制步长系数, [R1 R2]T为不平衡力 问题:外力不按比例加载?
其他数值方法
离散单元法 刚体弹簧元法 无网格法
其他数值方法
离散单元法 Cundall于1971年提出 基本原理:牛顿第二运动定律 建立每个刚体的力平衡方程, 采用动力松弛法求解 特点 岩体或颗粒组合体被模拟成通过边或角相互接触而产 生相互作用 块体之间边界的相互作用可以体现其不连续性和节理 的特征 使用显式积分迭代法,允许大位移、转动和非线性材料
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 割线迭代法
[δ 1 ] = [K 0 ]−1 [P]
[δ 2 ] = [K 1 ]−1 [P ] …… [δ n ] = [K n−1 ]−1 [P]
割线刚度 [K 1 ] 割线刚度 [K 2 ]
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 切线迭代法
有
δFs = K s ∆δ s
t + ∆t
Fs = t F + δFs ≤ c + t Fn tan ϕ
文献中给出了Ks和Kn的建议值
其他数值方法
运动方程 位移和转角、一阶导（速度）、二阶导（加速度） 采用中心差分法，逐步迭代，至达到平衡状态
考虑阻尼（Rayleigh阻尼）
阻尼系数的选择一般依赖经验方法
…
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 等刚度迭代法
[δ 1 ] = [K 0 ]−1 [P]
节点荷载 [P1 ] ，
[∆P1 ] = [P] − [P1 ]
[∆δ 2 ] = [K 0 ]−1 [∆P1 ]
[δ 2 ] = [δ 1 ] +其他数值方法
无网格法
无网格伽辽金法
由移动最小二乘法得到近似场函数 由变分原理得到基本方程
取
问题:适用范围?
刚度硬化 系数 0 ≤ λ ≤ 1
非线性方程组的解法
弧长法
[K ][∆u ] = ∆λ [P] + [R ]
i i +1 i +1 i
荷载控制 步长系数 同时控制 ∆λ 和 ∆u ,使得
不平衡力
[∆u ]T [∆u ] + ∆λ2 = dS
i +1
在 λ − u 空间中,向量r迭代按一圆弧进行,故称弧长法 ……可以求得 ∆λ
其他数值方法
时步选择 具有质量m和刚度系数k的弹性系统，其运动方程
时间步长需满足
∆t ≤ 2 m / k
块体系统由多个质量不一、刚度不同的块体组成， 理论上取其最小的 ∆tmin，实际计算时取小于 ∆tmin
其他数值方法
接触判断算法 二维相对成熟，三维有待发展 分格接触检索算法 1.将研究区域划分成网格 2.在每个块体内，按照接触算法，检索块体是否接触 3.每次迭代，均需接触检索
4. 区分不同情况,计算应力增量
非线性方程组的解法
弹塑性单元应力调整计算——二维
其他数值方法
其他数值方法
有限单元法的局限 基于连续体力学,相邻边界上位移协调, 对于处理应力或位移间断的问题麻烦. 目前解决办法: 1.重新划分网格 计算工作量大 2.改变单元位移插值函数使单元各向异性 网络依赖性,精度差 3.生死单元 网格太细密
[K n−1 ][∆δ n ] = [∆Pn ]
[∆δ n ] = [K n−1 ]−1 [∆Pn ]
[δ n ] = [δ n−1 ] + [∆δ n ]
步骤：1.施加第n步荷载增量 [∆Pn ]，利用始点刚度矩阵 [K n−1 ]，求位移增量 [∆δ n ] 2.计算单元应变增量和应力增量，计算总的位移与应力
非线性方程组的解法
弹塑性单元应力调整计算——一维
非线性方程组的解法
弹塑性单元应力调整计算——二维 1. 已知j-1步迭代结束后的所有量 2. 按弹性材料试算j步各单元高斯点处应力增量和 应力水平 3. 检查各单元高斯点处应力水平: j-1步弹性, j步弹性 j-1步屈服, j步(弹性)卸载 j-1步弹性, j步屈服 j-1步屈服, j步应力继续增长
其他数值方法
质点－桁架模型 弹簧刚度 轴向刚度等效
k1 = 1 Et 2(1 + ν )
ν k2 = Et 2 1 −ν
切向刚度等效
k2 = 1 Et 2(1 + ν )
取 ν = 1/ 3
开裂前，桁架单元的应力－应变关系服从混凝土 单轴应力应变关系
其他数值方法
两个例子
其他数值方法
两个例子
其他数值方法
非线性方程组的解法
非线性方程组的解法
主要有两类： 几何非线性，材料非线性 几何非线性 物理位置：平衡方程必须按照变形后的几何位置建立 数学描述：物体的变形一般由位移的一阶微分求得， 当变形很大，高阶微分不能忽略 举例：
非线性方程组的解法
主要有两类： 几何非线性，材料非线性 材料非线性 混凝土，钢筋，粘结滑移，时变效应 一般认为，混凝土开裂即存在材料非线性 其他非线性：边界非线性 常见的是接触问题
[K ][∆δ ] = [∆P ]
n n n
[∆δ n ] = [K n ]−1 [∆Pn ]
[δ n ] = [δ n−1 ] + [∆δ n ]
步骤：……. 缺点：存储两个刚度，存贮量大
非线性方程组的解法
逐步增量法－中点刚度增量法
施加一半荷载，求出相应的位移（中点位移）
[K n−1 ]∆δ n − 1 = 1 [∆Pn ]
非线性方程组的解法
逐步搜索法 可求出极限荷载，无法求解
步骤：1. 加一级荷载 ∆P ，结果发散，退回 2. 加 1 ∆P ，若结果发散，再次退回
2
P −δ
下降段
3. 加
1 ∆P ，……至收敛，得到极限荷载 4
非线性方程组的解法
虚加刚性弹簧法 加虚拟弹簧, 改变原结构形式，新的结构形式 曲线在一定范围内没有下降段
[B] = [M ][R][Q]
其他数值方法
刚体弹簧元
k n [D] = 0
Tresca屈服条件
0 ks
本构矩阵 (2x2)
平面应变下的莫尔-库仑屈服条件
平面应力下的莫尔-库仑屈服条件
其他数值方法
刚体弹簧元
其他数值方法
无网格法
特点 在计算域上，用离散的点拟合场函数，采用基于节点（非单元节点） 的位移插值函数 部分或彻底地消除网格，避免了有限单元法中单元网格的限制 在处理大变形、开裂等问题时，有其优势 采用离散的点拟合场函数的方法 核（权）函数法 移动最小二乘法 径向基函数法 然后建立方程求解
[δ 1 ] = [K 0 ]−1 [P]
节点荷载 [P1 ] ， 即时切线刚度 [K 1 ]
[∆P1 ] = [P] − [P1 ]
[∆δ 2 ] = [K1 ]−1 [∆P1 ]
[δ 2 ] = [δ1 ] + [∆δ 2 ]
[δ n ] = [δ n−1 ] + [∆δ n ]
节点荷载 [P2 ] ， 即时切线刚度 [K 2 ]
非线性方程组的解法
添加行列法
非线性方程组的解法
强制迭代法 结构总刚是负定的,只要施加负荷载,可求得相应的位移增量
问题:适用范围?
非线性方程组的解法
硬化刚度法 人为增大结构刚度,改善收敛性 一般结构切线刚度可分为
KT = K E − K p
弹性 非线性影响 刚度 刚度
K T = K E − λK p
能量准则
([P ] [δ ]) res
T
两曲线的接近如何描述？
非线性方程组的解法
两个问题： 考虑结构负刚度的一些算法 弹塑性单元应力调整计算
非线性方程组的解法
考虑结构负刚度的一些算法 问题：混凝土结构 P − δ 下降段，由力求位移求解困难 原因：上升段，刚度矩阵正定，线性方程组有解 下降段，刚度矩阵非正定，求解困难 负刚度：非正定的刚度矩阵 几种方法：逐步搜索法，虚加刚性弹簧法，位移加载法， 强制迭代法，硬化刚度法，弧长法
刚体弹簧元 刚体＋弹簧 自由度少，总刚阶数少 可在在极限分析领域运用 思路 弹簧单元刚度矩阵 界面应力应变关系 界面法向和切向的应力应变关系 界面法向和切向刚度
其他数值方法
刚体弹簧元
[K ] = ∫ [B ] [D ][B ]ds
T
k n [D] = 0
0 ks
单元刚度矩阵 (6x6) 本构矩阵 (2x2) 形函数矩阵 (2x6)
网格大小与数目的选择 块体大小，计算精度
其他数值方法
公共面法 块体间接触判断 公共面的定义（如何得到） 假象面 令两块体向薄板靠近 两块体在接触时停止运动 薄板被固定在一个位置，该面即为公共面 块体与公共面间的接触关系
其他数值方法
公共面法 接触类型
其他数值方法
颗粒离散元法 多边形、多面体 圆形、球形 可以以链杆或弹簧连接颗粒 质点－桁架模型
2
2
δ 1 = [δ n −1 ] + ∆δ 1 n− 2 n− 2
K n− 1 确定相应的刚度（中点刚度） 2
以该刚度作为初始刚度求当前位移增量
K n − 1 [∆δ n ] = [∆Pn ] 2
优点: 比两步欧拉折线法精度高
P −δ
问题:在何处加弹簧? 弹簧刚度如何?
非线性方程组的解法
位移控制法 支座位移控制加载,可求支座的反力, 和全部 P − δ 曲线 两个算法:刚度矩阵重新排列法, 添加行列法
非线性方程组的解法
添加行列法
[p1 p2]T为参考矢量, ∆λ 为荷载控制步长系数, [R1 R2]T为不平衡力 问题:外力不按比例加载?
其他数值方法
离散单元法 刚体弹簧元法 无网格法
其他数值方法
离散单元法 Cundall于1971年提出 基本原理:牛顿第二运动定律 建立每个刚体的力平衡方程, 采用动力松弛法求解 特点 岩体或颗粒组合体被模拟成通过边或角相互接触而产 生相互作用 块体之间边界的相互作用可以体现其不连续性和节理 的特征 使用显式积分迭代法,允许大位移、转动和非线性材料
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 割线迭代法
[δ 1 ] = [K 0 ]−1 [P]
[δ 2 ] = [K 1 ]−1 [P ] …… [δ n ] = [K n−1 ]−1 [P]
割线刚度 [K 1 ] 割线刚度 [K 2 ]
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 切线迭代法
有
δFs = K s ∆δ s
t + ∆t
Fs = t F + δFs ≤ c + t Fn tan ϕ
文献中给出了Ks和Kn的建议值
其他数值方法
运动方程 位移和转角、一阶导（速度）、二阶导（加速度） 采用中心差分法，逐步迭代，至达到平衡状态
考虑阻尼（Rayleigh阻尼）
阻尼系数的选择一般依赖经验方法
…
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 等刚度迭代法
[δ 1 ] = [K 0 ]−1 [P]
节点荷载 [P1 ] ，
[∆P1 ] = [P] − [P1 ]
[∆δ 2 ] = [K 0 ]−1 [∆P1 ]
[δ 2 ] = [δ 1 ] +其他数值方法
无网格法
无网格伽辽金法
由移动最小二乘法得到近似场函数 由变分原理得到基本方程
取
问题:适用范围?
刚度硬化 系数 0 ≤ λ ≤ 1
非线性方程组的解法
弧长法
[K ][∆u ] = ∆λ [P] + [R ]
i i +1 i +1 i
荷载控制 步长系数 同时控制 ∆λ 和 ∆u ,使得
不平衡力
[∆u ]T [∆u ] + ∆λ2 = dS
i +1
在 λ − u 空间中,向量r迭代按一圆弧进行,故称弧长法 ……可以求得 ∆λ
其他数值方法
时步选择 具有质量m和刚度系数k的弹性系统，其运动方程
时间步长需满足
∆t ≤ 2 m / k
块体系统由多个质量不一、刚度不同的块体组成， 理论上取其最小的 ∆tmin，实际计算时取小于 ∆tmin
其他数值方法
接触判断算法 二维相对成熟，三维有待发展 分格接触检索算法 1.将研究区域划分成网格 2.在每个块体内，按照接触算法，检索块体是否接触 3.每次迭代，均需接触检索
4. 区分不同情况,计算应力增量
非线性方程组的解法
弹塑性单元应力调整计算——二维
其他数值方法
其他数值方法
有限单元法的局限 基于连续体力学,相邻边界上位移协调, 对于处理应力或位移间断的问题麻烦. 目前解决办法: 1.重新划分网格 计算工作量大 2.改变单元位移插值函数使单元各向异性 网络依赖性,精度差 3.生死单元 网格太细密
[K n−1 ][∆δ n ] = [∆Pn ]
[∆δ n ] = [K n−1 ]−1 [∆Pn ]
[δ n ] = [δ n−1 ] + [∆δ n ]
步骤：1.施加第n步荷载增量 [∆Pn ]，利用始点刚度矩阵 [K n−1 ]，求位移增量 [∆δ n ] 2.计算单元应变增量和应力增量，计算总的位移与应力
非线性方程组的解法
弹塑性单元应力调整计算——一维
非线性方程组的解法
弹塑性单元应力调整计算——二维 1. 已知j-1步迭代结束后的所有量 2. 按弹性材料试算j步各单元高斯点处应力增量和 应力水平 3. 检查各单元高斯点处应力水平: j-1步弹性, j步弹性 j-1步屈服, j步(弹性)卸载 j-1步弹性, j步屈服 j-1步屈服, j步应力继续增长
其他数值方法
质点－桁架模型 弹簧刚度 轴向刚度等效
k1 = 1 Et 2(1 + ν )
ν k2 = Et 2 1 −ν
切向刚度等效
k2 = 1 Et 2(1 + ν )
取 ν = 1/ 3
开裂前，桁架单元的应力－应变关系服从混凝土 单轴应力应变关系
其他数值方法
两个例子
其他数值方法
两个例子
其他数值方法
非线性方程组的解法