非线性动力学方程的求解方法
非线性动力系统的数值计算方法及稳定性分析
非线性动力系统的数值计算方法及稳定性分析非线性动力系统是一种非常常见的实际物理系统，例如电路、化学反应、天气系统等，它们的行为通常比线性系统更加复杂。
数值计算非线性动力系统的稳定性与动力学特性是一个非常重要的课题，对于研究和预测实际系统的行为有着非常重要的意义。
在本文中，我们将介绍几种常见的非线性动力系统的数值计算方法及它们的稳定性分析。
一、欧拉法欧拉法是动力系统数值计算中最基本的一种方法。
它的基本思路是将连续的时间离散化，将微分方程转化成差分方程，然后用迭代的方式求解。
欧拉法的迭代公式为：$$y_{n+1}=y_{n}+hf(y_n)$$其中，$h$为步长，$f(y_n)$是微分方程在$y_n$处的导数。
欧拉法是一种比较简单易懂的方法，但是它的稳定性较差，容易产生数值误差。
欧拉法对于初始值的依赖性很强，如果步长$h$选取过大，就会导致解的不稳定。
因此，在使用欧拉法进行数值计算时，我们需要根据实际问题来调整步长，以保证数值解的正确性。
二、龙格-库塔法龙格-库塔法是一种常见的数值积分方法，在动力系统数值计算中也常常被使用。
它的基本思路是利用微分方程的某些性质，选取合适的时间步长和权重，在数值上求得微分方程的积分近似解。
龙格-库塔法通常可以由一些权重系数和步长系数组成，如下：$$Y_{n+1}=Y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$$$$k_1=hf(Y_n)$$$$k_2=hf(Y_n+\frac{1}{2}k_1)$$$$k_3=hf(Y_n+\frac{1}{2}k_2)$$$$k_4=hf(Y_n+k_3)$$其中，$k_1,k_2,k_3,k_4$均为微分方程在相应位置处的导数。
龙格-库塔法比欧拉法更加稳定，适用于多数动力系统的数值计算。
但是，龙格-库塔法在计算一些比较长时间范围内的运动时，需要降低步长以保证解的精度。
同时，权重系数和步长系数也需要根据具体问题调整，才能得到更加准确的数值解。
解非线性方程的牛顿迭代法及其应用
解非线性方程的牛顿迭代法及其应用一、本文概述非线性方程是数学领域中的一个重要研究对象，其在实际应用中广泛存在，如物理学、工程学、经济学等领域。
求解非线性方程是一个具有挑战性的问题，因为这类方程往往没有简单的解析解，需要通过数值方法进行求解。
牛顿迭代法作为一种古老而有效的数值求解方法，对于求解非线性方程具有重要的应用价值。
本文旨在介绍牛顿迭代法的基本原理、实现步骤以及在实际问题中的应用。
我们将详细阐述牛顿迭代法的基本思想，包括其历史背景、数学原理以及收敛性分析。
我们将通过具体实例，展示牛顿迭代法的计算步骤和实际操作过程，以便读者能够更好地理解和掌握该方法。
我们将探讨牛顿迭代法在各个领域中的实际应用，包括其在物理学、工程学、经济学等领域中的典型应用案例，以及在实际应用中可能遇到的问题和解决方法。
通过本文的介绍，读者可以深入了解牛顿迭代法的基本原理和应用技巧，掌握其在求解非线性方程中的实际应用方法，为进一步的研究和应用提供有力支持。
二、牛顿迭代法的基本原理牛顿迭代法，又称为牛顿-拉夫森方法，是一种在实数或复数域上近似求解方程的方法。
其基本原理是利用泰勒级数的前几项来寻找方程的根。
如果函数f(x)在x0点的导数f'(x0)不为零，那么函数f(x)在x0点附近可以用一阶泰勒级数来近似表示，即：这就是牛顿迭代法的基本迭代公式。
给定一个初始值x0，我们可以通过不断迭代这个公式来逼近f(x)的根。
每次迭代，我们都用当前的近似值x0来更新x0，即：这个过程一直持续到满足某个停止条件，例如迭代次数达到预设的上限，或者连续两次迭代的结果之间的差小于某个预设的阈值。
牛顿迭代法的收敛速度通常比线性搜索方法快，因为它利用了函数的导数信息。
然而，这种方法也有其局限性。
它要求函数在其迭代点处可导，且导数不为零。
牛顿迭代法可能不收敛，如果初始点选择不当，或者函数有多个根，或者根是重根。
因此，在使用牛顿迭代法时，需要谨慎选择初始点，并对迭代过程进行适当的监控和调整。
牛顿迭代法解动力学方程不收敛
牛顿迭代法解动力学方程不收敛摘要：1.引言2.牛顿迭代法简介3.动力学方程及其收敛性问题4.牛顿迭代法在解动力学方程中的应用5.牛顿迭代法在解动力学方程中的不收敛问题6.结论正文：1.引言在物理学中，动力学方程是描述物体运动状态的数学模型，广泛应用于各种实际问题中。
然而，在求解动力学方程时，常常会遇到收敛性问题。
牛顿迭代法作为一种求解非线性方程的数值方法，被广泛应用于解动力学方程。
本文将探讨牛顿迭代法在解动力学方程中的不收敛问题。
2.牛顿迭代法简介牛顿迭代法是一种求解非线性方程的数值方法，其基本思想是通过迭代使得函数值逐步逼近零。
对于非线性方程F(x) = 0，牛顿迭代法的迭代公式为：x[n+1] = x[n] - F(x[n])/F"(x[n])，其中F"(x) 表示F(x) 的导数。
牛顿迭代法具有二阶收敛性，即当迭代步长足够小，且初始值足够接近真实解时，可以通过有限次迭代得到精确解。
3.动力学方程及其收敛性问题动力学方程描述了物体在给定力的作用下的运动状态，通常包括质量、速度、加速度等物理量。
求解动力学方程时，通常需要采用数值方法，因为解析解往往难以求得。
然而，在数值求解过程中，可能会遇到收敛性问题。
例如，在迭代过程中，如果迭代步长过大或者初始值与真实解差距过大，可能导致迭代结果发散，无法得到精确解。
4.牛顿迭代法在解动力学方程中的应用由于牛顿迭代法具有二阶收敛性，因此在求解动力学方程时，可以得到较好的数值解。
在实际应用中，可以根据动力学方程的特点，选择合适的牛顿迭代法求解。
例如，对于具有显式解的动力学方程，可以直接使用牛顿迭代法求解；对于具有隐式解的动力学方程，可以通过拟合等方法得到显式解，然后使用牛顿迭代法求解。
5.牛顿迭代法在解动力学方程中的不收敛问题尽管牛顿迭代法具有二阶收敛性，但在求解动力学方程时，仍然可能出现不收敛的情况。
这主要是因为动力学方程的非线性特性和迭代过程中的误差累积。
数值分析中的非线性方程求解与优化
数值分析中的非线性方程求解与优化在数值分析领域中，非线性方程求解是一个重要的问题。
许多实际问题都可以被建模为非线性方程，而求解这些方程对于解决实际问题具有重要意义。
本文将介绍非线性方程求解的基本概念、方法和优化技术。
一、非线性方程求解的概念非线性方程是指方程中包含非线性项的方程。
与线性方程不同，非线性方程的解不再是一条直线，而是一条曲线或曲面。
非线性方程的求解是寻找方程中满足特定条件的变量值或函数的过程。
二、非线性方程求解的方法1. 迭代法迭代法是解决非线性方程求解问题中常用的方法。
迭代法的基本思想是通过不断逼近方程的解，使得迭代序列逐步收敛于方程的解。
常见的迭代法包括牛顿迭代法、割线法和弦截法等。
以牛顿迭代法为例，假设要求解方程f(x) = 0，首先选择一个初始估计值x0，然后通过迭代公式进行迭代计算直到满足收敛条件。
迭代公式为：xn+1 = xn - f(xn)/f'(xn)，其中f'(xn)表示f(x)在xn处的导数。
2. 区间划分法区间划分法是通过将求解区间划分为若干个子区间，然后在每个子区间内搜索方程的解。
这种方法常用于求解具有多个解的非线性方程。
一般可以使用二分法、割线法和弦截法等算法进行区间划分和求解。
3. 优化技术优化技术常用于求解非线性方程的最优解。
在数值分析中，优化问题可以理解为寻找使得目标函数达到最大或最小值的变量值。
常用的优化算法包括梯度下降法、拟牛顿法和粒子群算法等。
这些算法通过迭代过程不断调整变量值，使得目标函数逐渐趋于最优解。
三、非线性方程求解与优化的应用非线性方程求解和优化技术在实际问题中具有广泛的应用。
以下是一些应用领域的例子：1. 工程领域：在工程设计中，需要求解非线性方程以确定优化的设计参数。
例如，在机械设计中，可以通过求解非线性方程来确定零件的几何尺寸和运动轨迹。
2. 金融领域：在金融衍生品定价和风险管理中，需要求解非线性方程来估计资产价格和风险敞口。
非线性动力学
t∈R
x∈ Rn
的解，则显然它是不仅是时间的函数，而且也是初值的函数，即解随着初值的改变而改变， 可以将解记为
φ(t, x0 )
当 x0 是 R n 中的某一点时，φ (t, x0 ) 代表了 1 条解轨线，而
{φ(t, x0 ) x0 ∈ D}
则代表了一族轨线。将φ看成是一个映射，即
φ : R× Rn → Rn
运动行为，它在物理上对应了这样的一个观点：在系统的最初阶段，系统由于外界的初始干 扰，将呈现相当复杂的运动形式，但随着时间的延续，运动将进入平稳状态，而这种平稳状 态体现了动态系统的本质结构。
微分方程解的最终形态通常有： (1) 平衡点 (2) 周期解 (3) 拟周期解 (4) 混沌解
6.4.1 平衡点
图 6-7 所示是 2 维线性系统的相轨线，坐标原点是系统的平衡点，图 6-7a、b 中的平衡 点是稳定的，称为稳定结点，图 6-7c 中的平衡点是不稳定的，称为鞍点。
图 6-7 2 维线性系统的相轨线
6.5.2 任意解的稳定性
设 x = ψ (t)是微分方程 x& = F(t, x)
第 6 章 非线性动力学
-0.5
-1
-1.5
0.5
1
1.5
图 6-2 例 1 相图
例2
如图 6-3 所示是微分方程
&y& + 0.2 y& + y = 0
在相平面 (x1, x2 ) ，
x1 = y
x2 = y&
上的轨线图，平衡点为 (0,0)，当 t → ∞ 时，解轨线趋于平衡点。
0.6 0.4 0.2
-0.6
-0.4
-0.2 -0.2
使用matlab求解vanderpol方程的研究方法
使用matlab求解vanderpol方程的研究方法研究van der pol方程是一种重要的非线性动力学问题，它描述了一些自振系统的运动行为。
这个方程由物理学家布里尔·范·德·波尔在1920年首次提出，并广泛用于描述电路、生物和化学系统中的自振现象。
van der pol方程的一般形式如下：d^2x/dt^2 - μ(1 - x^2)dx/dt + x = 0其中，x是系统的状态变量，t是时间，而μ是一个非负常数，代表了系统的非线性因素。
为了研究van der pol方程的行为，我们可以使用MATLAB进行求解和分析。
下面将介绍几种常用的方法。
1. 数值解法：通过数值方法求解van der pol方程是最常用的研究方法之一、MATLAB中可以使用ode45函数来进行求解。
该函数基于Adams-Bashforth-Moulton方法，可以自动选择适当的步长，并给出较高的数值精度。
例如，我们可以定义一个匿名函数来表示van der pol方程，并使用ode45来求解：```MATLABmu = 1; % 设置μ的值[t, x] = ode45(f, [0, 10], [0.1; 0]); % 求解方程plot(t, x(:,1)); % 绘制解的图像xlabel('t'); ylabel('x');```这段代码定义了一个匿名函数f，其中t是时间，x是状态变量。
通过ode45函数对该方程进行求解，给定了初始条件[0.1; 0]，求解时间范围为[0, 10]。
最后，使用plot函数对解进行可视化。
2. 分岔和稳定性分析：van der pol方程是一个具有非线性耗散的系统，因此它可以表现出各种动力学行为，如周期振荡、混沌和吸引子等。
MATLAB提供了一些函数用于分析方程的分岔和稳定性。
例如，可以使用bifurcation函数分析van der pol方程的分岔图：```MATLABmu = linspace(0, 10, 100); % 设定μ的取值范围x = linspace(-2, 2, 100); % 设定x的取值范围bifplot(bif_x, bif_par); % 绘制分岔图xlabel('μ'); ylabel('x');```这段代码根据设定的μ和x范围，使用bifurcation函数计算系统的分岔图，并使用bifplot函数进行可视化。
第十一章 非线性动力学
可饱和的代谢过程；酶诱导；较高剂量时 的肝中毒；肝血流的变化；代谢物的抑制 作用
二、非线性药物动力学特点与识别
特点：
药物消除为非一级动力学，遵从米氏方程 AUC与剂量不成正比 消除半衰期随剂量增大而延长，剂量增加至一定 程度时，半衰期急剧增大 动力学过程可能会受到合并用药的影响 代谢物的组成比例受剂量的影响
当C0>>Km时, t1/2=C0/(2Vm) 当Km>>C0时, t1/2=0.693Km/Vm
清除率Cl
dX dt Cl C VmC dX dt ( dC dt ) V V Km C Vm V Cl Km C
当C>>Km时, Cl与C成反比：CL=Vm*V/C 当Km>>C时, Cl与C无关： CL=Vm*V/Km
线性动力学
血药浓度与剂量呈正比 ； AUC与剂量呈正比；t1/2、k、 V、Cl与剂量无关
非线性动力学
Dose-dependant PK 动力学参数与剂量有关 存在饱和现象
k
AUC
t1/2
X0
X0
X0
注：图中实线表示非线性，虚线表示线性非线性药代动力学主要见于：
与药物代谢有关的可饱和的酶代谢过程； 与药物吸收、排泄有关的可饱和的载体转 运过程； 与药物分布有关的可饱和的血浆/组织蛋白 结合过程； 酶诱导及代谢产物抑制等其他特殊过程。
五、非线性动力学参数的求算
1. Km及Vm的求算：根据-dC/dt 求算
dC Vm C dt K m C
Lineweaver-Burk方程式： Hanes-Woolf方程式： Eadie-Hofstee方程式：
求解非线性动力学方程的预测_校正数值算法
求解非线性动力学方程的预测2校正数值算法黄涛樊建平何建平王乘( 华中科技大学工程计算与仿真研究所,武汉,430074)摘要提出了一种求解非线性动力学方程的预测2校正数值算法.由于校正过程消除了线性化带来的误差,使计算精度得到较大提高.数值算例表明该方法正确、有效.关键词非线性动力学, 直接积分法, 内共振, 饱和(渗透) 现象, 跳跃现象, 分岔0 引言求解非线性动力学方程,有解析法和数值法. 虽然解析法在不断发展, 但迄今为止, 大多数非线性动力学方程仍然难以得到解析解, 为了满足科学研究和工程实际需要,必须借助数值方法.传统的非线性动力学问题的数值解法是直接积分法, 如中心差分法、Wil s o n2θ法、Ho u bolt 法和New ma r k 法等. 这类算法通过构造各种差分格式, 推导出前一时刻和后一时刻的递推关系, 再由初始条件求出各相继离散时刻的解. 文献[1 ]称此类算法为差分型直接积分法. 此类算法虽然得到广泛应用, 但其计算精度不高一直是个突出问题.另一类直接积分法为积分型直接积分法[1 ] , 该方法在时间步长上,对状态空间中的非线性动力学方程进行线性化, 将问题转化为求解线性动力学方程问题,其实质是用线性方程的解近似非线性方程的解. 影响此类算法计算精度的主要原因: (1) 非线性方程线性化的误差; (2) 线性方程本身的计算误差. 在钟万勰提出指数矩阵的精细计算法[2 ,3 ] 后,对于线性动力学方程´v= B v + g ( t) ( 1) 式中B 是非奇异常数矩阵, g( t) 是t 的多项式,可以求得在计算机有效字长范围内的精确解. 因此对于非线性方程的积分型直接积分法, 只有线性化的精度是影响其计算精度的主要原因.文献[4 ]2[6]分别提出了不同的线性化方法. 通过比较这三种不同的方法可以看出:文献[4 ] 、[5 ]提出的方法, 矩阵 B 只在分段区间上是常矩阵, 不同的分段区间B 的值不同,这样不仅矩阵有可能出现奇异,而且增加了计算量,降低了计算效率. 文献[6 ]的方法, B 矩阵在整个求解区间都是不变的常矩阵,避免了上述问题.文献[4 ]2[6 ]还提出了提高计算精度的方法,文献[ 5 ]通过修正矩阵B 和一次多项式g ( t) 的系数来提高线性化的精度,从而提高非线性方程求解的精度; 文献[ 6 ]通过提高多项式g ( t) 的阶次以更准确逼近非线性方程, 从而达到提高计算精度的目的. 文献[6 ]的方法比文献[5]的方法计算精度高,但为了达到高精度,文献[6]的方法需要计算高阶导数.文献[ 5 ] 、[ 6 ]的方法都不能从根本上消除用线性方程代替非线性方程所产生的离散误差. 为了消除这一误差, 文献[ 4] 提出了一种不同的思路, 即仅把( 1) 式的解作为v 的预测值,通过重构非线性系统等价方程, 采用Newto n2Rap h so n 迭代求得v 的校正值, 由于此方法本身不存在离散误差, 从而使计算精度得到较大提高.本文基于文献[ 4]的思路, 构造了v 的三次插值多项式, 并结合文献[ 6 ] 提出的线性化方法, 提出了一种新的高精度非线性动力学方程预测2校正数值算法. 算例表明本文方法具有较高的计算精度.1 非线性动力学方程的预测解状态空间中n 维自治或非自治非线性系统可表示为´v= F(v , t) ( 2)将F( v , t) 分离成线性和非线性两部分,即´v= H v + f ( v , t) ( 3)+ +则由 i 次计算 i + 1 次指数矩阵变为矩阵乘积系数 的加法运算步长τ= 0 . 1 s , 指数矩阵截断阶次 r = 6 , 小区间数 m = 10000 , 用本文方法计算 v 1 , 并与文献 [ 6 ] 、[ 4 ] 所 eH ( i + 1)Δh= I + ( a + 1) ( H Δh ) + a + b + 1 2( H Δh ) 2+ 述方法和四阶 Runge 2 K ut t a 法进行比较. 文献 [ 6 ]中 的两种方法分别是 f ( v , t) 的近似多项式取一次式 a + b + c + 1 ( H Δh ) 3 + a b + c + d + 1 ·( j = 1) 和取三次式 ( j = 3) 的情形 ; 文献 [ 4 ]的方法取 2 6 6 2 24 本文相同的参数 r = 6 , m = 10000 ; 四阶 Runge 2 K ut 2( H Δh ) 4 + a + b + c + d + e + 1 ( H Δh ) 5 +t a 的步长分别取 0 . 1 和 0 . 001 . 计算结果如表 1 所24 6 2 120示 .a+ b + c + d + e + f +1( H Δh ) 6 ( 15)120 24 6 2 720 以四阶 Runge 2 K ut t a 法 , 步长 0 . 001 计算的 v 1 此算法为无条件稳定算法.至此 , 对于 ( 3) 式表示的非线性动力学方程 , 由 ( 8) 式求其预测值 , 再由 ( 12) 式求其校正值 , 经过校 正后的解有较高的计算精度.3 算例例 1 单摆的非线性振动 . 运动微分方程为θ¨ ω2 si n θ= 0 , 做变换 v 1 =θ, v 2 = ´v 1 , 相应的一阶微分方 程组 ´v = H v + f ( v , t ) 为 作参考 , 可以看出本文方法计算的 v 1 最 接 近精 确 解 . 文献 [ 6 ]中 j = 1 时算法的解 , 相当于本文没有经 过校正的预测解 , 从结果可以看出由于线性化的误 差 , 求解时间的越长 , 解的误差越大 , 但是经过本文 方法校正后 , 不仅解的精度有较大提高 , 而且随着时 间的增加 , 解依然保持较高的精度 , 从而证明了本文 方法的有效性 . 与文献 [ 6 ]中 j = 3 的算法相比 , 本文方法计算的 v 1 精度还要稍高 , 但却避免了求 f ( v , t ) 的高阶导数 , 因而具有一定的优越性. 与文献 [ 4 ] 方 法相比 , 本文方法不仅精度高 , 而且由于 B 矩阵在´v 1=´v 2 0 1 - ω2 0 v 1+v 2ω2( v 1 - si n v 1 )( 16)整个求解区间都是常数矩阵 , 避免了重复计算指数矩阵 , 因而计算效率也高.取 ω2= 9 . 80665 , 初始条件 v 1 = 1 . 04720 , v 2 = 0 . 0 ,表 1 各种方法计算的单摆非线性振动幅角的比较时间/ s 本文方法( 步长 0 . 1) 文献 [ 6 ]方法( j = 1)( 步长 0 . 1) 文献 [ 6 ]方法( j = 3)( 步长 0 . 1) 文献 [ 4 ]方法( 步长 0 . 1) 四阶Run ge 2 Kut t a ( 步长 0 . 1) 四阶Run ge 2 Kut t a ( 步长 0 . 001) 1 - 1 . 022 33 - 1 . 021 41 - 1 . 022 33 - 1 . 021 75 - 1 . 022 24 - 1 . 022 33 2 0 . 948 47 0 . 945 51 0 . 948 46 0 . 947 54 0 . 948 25 0 . 948 47 3 - 0 . 828 00 - 0 . 822 01 - 0 . 827 99 - 0 . 827 16 - 0 . 827 62 - 0 . 827 99 4 0 . 665 34 0 . 655 52 0 . 665 34 0 . 665 20 0 . 664 81 0 . 665 31 5 - 0 . 467 36 - 0 . 453 35 - 0 . 467 40 - 0 . 468 64 - 0 . 466 72 - 0 . 467 33 6 0 . 243 66 0 . 225 73 0 . 243 75 0 . 247 03 0 . 242 98 0 . 243 61 7 - 0 . 006 19 0 . 014 58 - 0 . 006 32 - 0 . 012 17 - 0 . 005 55 - 0 . 006 14 8 - 0 . 231 63 - 0 . 253 49 - 0 . 231 46 - 0 . 222 83 - 0 . 232 14 - 0 . 231 69 9 0 . 456 36 0 . 477 17 0 . 456 18 0 . 444 95 0 . 456 67 0 . 456 42 10 - 0 . 655 92 - 0 . 673 63 - 0 . 655 76 - 0 . 642 54 - 0 . 656 00 - 0 . 655 97 20 - 0 . 255 64 - 0 . 222 51 - 0 . 255 97 - 0 . 300 51 - 0 . 256 49 - 0 . 255 50 50 0 . 195 37 0 . 133 51 0 . 194 15 - 0 . 108 78 0 . 175 33 0 . 195 71 100 - 0 . 980 68 - 0 . 991 88 - 0 . 981 72 - 0 . 734 32 - 1 . 004 14 - 0 . 980 44 5000 . 276 920 . 516 43 0 . 216 52 0 . 243 63- 0 . 790 170 . 280 26F124ω23 例 2 弹簧摆的受迫振动 . 如图 1 所示 , 设弹簧 摆存在微弱的阻尼 , 轴向阻尼力与质点的轴向速度 成正比 , 切向阻尼力与弹簧摆动角速度成正比 , 摆上 沿轴向作 用 一 个 简 谐 激 励 力 . 系 统 的 运 动 微 分 方 程[ 7 ] 为¨x + 2εc 1 ´x + kx + g ( 1 - co s θ) - ( l + x )θ´2 =mε2m co s (Ωt )θ¨+ 2εc 2θ´+g l + xsi n θ+2 l + x´x θ´= 0( 17)令 ω2 = k/ m ,ω2= g/ l , 做变换 v 1 = x , v 2 =θ, v 3 = ´v 1 ,1 2 v 4 = ´v 2 , 相应的一阶微分方程组 ´v = H v + f ( v , t) 为图 2 v 1, v 2时间历程曲线 (F = 6)´v 1´v 2 =´v 3 ´v 40 0 1 0 v 10 01 v 2+- ω2- 2εc 1v 0- ω20 - 2εc 2v 40 0ε2 F- g ( 1 - co s v 2 ) + ( l + v 1 ) v 2+co s (Ωt )m2 v 2 -gl + v 1 si n v 2 -2l + v 1v 3 v 4( 18)图 1 受轴向激励的带阻尼的弹簧摆设 k = 49 、m = 1 、g = 9 . 8 、l = 0 . 8 、Ω = 7 . 1 、c 1 =1 、c2 = 1 和ε= 0 . 05 , 则 ω1 = k/ m = 7 ,ω2 =g/ l = 3 . 5 , 由于 ω1 = 2ω2 ,由非线性振动理论可知 , 该弹簧 摆系统存在内共振. 在初值 v 1 = 0 , v 2 = 5°, v 3 = 0 , v 4= 0 下 , 取步长τ= 0 . 02 s , 其余参数与例 1 相同 , 分别用激励振幅 F = 6 , F = 12 , F = 13 , F = 20 , 分析该 系统的运动 , 计算结果示于图 2 2图 9 中 .图 3 v 1 , v 2 时间历程曲线 ( F = 12)图 4 v 1 , v 2 时间历程曲线 ( F = 13)比较 F = 6 、F = 12 两种情况可以看出 , 弹簧摆 轴向运动的稳态响应幅值随着激励 F 的增大而增 大 , 且和 F 近似呈线性关系 ; 而由于阻尼的存在 , 弹 簧摆切向运动最终趋于稳定焦点 ( 0 , 0) .图 5 v 1 , v 2 时间历程曲线 ( F = 20)图 8 弹簧摆切向运动相平面图 ( F = 12)图 6 弹簧摆轴向运动相平面图图 7 弹簧摆切向运动相平面图 ( F = 6)比较 F = 13 、F = 20 两种情况可以看出 , 弹簧摆 轴向运动的稳态响应幅值不再 随着 激励 增大 而 增 大 , 而是保持不变 , 出现了“饱和现象”; 相反 , 弹簧摆 切向运动不再趋于焦点 ( 0 , 0) , 而是出现了稳态响 应 , 且幅值随着 F 增大而增大 , 表明轴向激励输入 的能量向切向自由度发生了转 移 , 出 现 了“渗 透 现 象”.图 9 弹簧摆切向运动相平面图比较 F = 12 、F = 13 两种情况可以看出 , 在弹簧 摆轴向运动幅值变化不大的情况下 , 弹簧摆切向运 动稳态幅值突然从 0 变化到一个有限值 , 发生了“跳 跃现象”. 进一步的计算分析表明 , 在区间 ( 12 . 7 , 12 . 8) 内 , 存在参数 F 的一个动态分岔点 , 当激励 F 不断 增大经过此分岔点时 , 就会发生跳跃现象和饱和 ( 渗 透) 现象 . 在步长τ= 0 . 02 s 的情况下 , 用本文方法 计算得出的结果与理论分析的结论[ 7 , 8 ] 一致 , 说明 了本文方法的有效性和可靠性.4 结论本文提出了一种求解非线性动力学方程的预测 2校正数值算法 , 该方法线性化过程简单 , 由于对 预测解的校正消除了线性化带来的离散误差 , 使得 计算精度得到较大提高 , 算例表明 , 即使采用较大的时间步长(如τ= 0 .1 s) 和较长的求解时间(如t = 500 s) , 用本文方法求得的解依然保持较高的精度.整个计算过程用矩阵表示, 算法简单, 易于编写程序, 可用于多自由度、强非线性, 非保守系统.参考文献[ 1 ]吕和祥, 于洪洁, 裘春航.动力学方程的积分型直接积分法.应用数学和力学, 2001 , 22 (2) : 1512156 . [ 2 ] 钟万勰, 杨再石.连续时间L Q 控制主要本征对的算法.应用数学和力学, 1991 , 12 (1) : 45250 .[ 3 ]钟万勰.暂态历程的精细计算方法.计算结构力学及其应用, 1995 , 12 (1) : 125 .[ 4 ]郑兆昌, 沈松, 苏志霄.非线性动力学常微分方程组高精度数值积分方法.力学学报, 2003 , 35 (3) : 2842295 .[ 5 ]张洵安, 姜节胜.结构非线性动力方程的精细积分算法.应用力学学报, 2000 , 17 (4) : 1642168 .[ 6 ]裘春航, 吕和祥, 钟万勰.求解非线性动力学方程的分段直接积分法.力学学报, 2002 , 34 (3) : 3692378 . [ 7 ]刘延柱, 陈立群.非线性振动.北京: 高等教育出版社, 2001 .[ 8 ] 奈弗A H ,穆克D T著.非线性振动.宋家骕, 罗惟德, 陈守吉译.北京: 高等教育出版社, 1990 .A P R E D ICT2CO RRECT NUM ERICA L I NTEGRATIO N SC HEM E FO RSOL VI NG NO N L I N EA R DY NAM IC EQUATIO N SH u a n g Tao Fa n J ia n pi n g He J ia n pi n g Wa n g Che n g( En g i nee r i n g Com p ut a t i o n a n d S i m u l at i o n I n st i t u d e , H u a z hon g V ni v e r s i t yo f S ci e nce a n d Tec h onol o g y , W u h a n , 430074)Abstract A new nu me rical i nt egratio n sche me i nco rpo rati ng a p redict2co r rect al go rit h m fo r sol v i n g t h e no nli nea r dyna mic sy st e ms wa s p ropo sed i n t h i s p ap e r . Si nce t he e r ro r ca u sed by li nea rizati o n ca n be e2 li mi nat e d i n t he co r rectio n p roce ss , t h e calculatio n accuracy i s great l y i mp ro ve d. N u mericl re s ult s sho w t h e co r r ect n e s s a n d vali d it y of t h e met h o d.K ey w ords no n li n ea r dyna m ic s , di r ect i n t e g ratio n met h o d , i n t e r n al re s o n a n ce , sat u ratio n p h e n o m e2 no n , j u mpi n g p h e n o me n o n , bif u rcatio n。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
非线性动力学方程的求解方法1、概述在工程实际问题中，我们常常面临这样的选择：我们所遇的问题究竟是静力的还是动力的。
静力问题与动力问题，从力学的角度看就是是否考虑与加速度有关的力，而从数学求解方法看则是一个三维边值问题还是一个四维边值-初值问题。
在这个问题的选择上没有固定的原则，一般取决于我们研究者、分析者对工程问题的判断。
一般认为，实际工程大都是处于动力环境之中，因而属于动力问题。
但是，由于时间、经费等方面的原因的限制，我们不可能把所有的问题都按照动力问题的方法来分析。
对于许多具体的问题，与速度和加速度有关的力足够小，但是又影响结构分析结果的，将采用静力假定来模拟这些力。
线性的动力有限元控制方程如式(1-1)所示。
[]}{}]{[}]{[}{R q K q D qM =++ (1-1) 式中[M ][D ][K ]分别为结构的质量、阻尼和刚度矩阵，{R }为荷载列矢量，}{q、}{q 和}{q 分别是加速度、速度和位移列矢量。
式(1-1)的解法大体上可以分为两类：直接积分法和模态叠加法。
直接积分法在对控制方程进行数值积分之前不对方程做任何形式的变换，直接用数值积分的方法在时域上一步一步地对方程进行积分。
模态叠加法是在求解之前对方程进行某种数学变换，使基底降低，或使矩阵的带宽减小，再进行求解。
这两种方法在形式上不同，但是密切相关。
上述每一类求解方法中又有许多具体的解法，每一种解法又有各自的特点。
因此我们在选择一种方法求解一个问题时，要对该方法的收敛性、稳定性、效率、精度和费用等进行一些分析，讨论它对所求问题的有效性，从而使我们能够针对某一特定的问题，选择合适的方法。
直接积分法基于以下两条：(1)不是在求解时间区间内任意时刻t 都满足式(1-1)，而是在相隔△t 上的一些离散时刻满足式(1-1)。
(2)对位移、速度和加速度在每一时间区间△t 内变化的形式进行假设，事实上若把式(1-1)看成一个常系数微分方程组，便可以用任何一种有限差分格式通过位移来近似表示速度和加速度，因此不同的差分格式就得到不同的方法。
从差分格式上看，分显式与隐式两大类方法，所谓显式差分法就是不必对方程进行求解，而是由前一时刻t 的平衡条件直接可求解△t 增量后t+△t 时刻的各参数解。
而隐式差分法则必须对方程进行求解。
显式差分法有中心差分法等。
隐式差分法有Wilson -θ法和Newmark 法。
不同方法解的精度、稳定性、收敛性、效率不同。
2、中心差分法2.1中心差分法简述中心差分法的差分格式为：}]{}{2}[{1}{2t t t t t t q q q tq ∆+∆-+-∆= (2-1) }]{}{[21}{t t t t t q q tq ∆+∆-+-∆= (2-2) 认为}{t q和}{t q 在任意时刻t 上满足平衡方程(1-1)，即： []}{}]{[}]{[}{t t t t R q K q D qM =++ (2-3) 将式(2-1)和式(2-2)代入到式(2-3)中，可得：}]){[21][1(}]){[1]([}{}]{[21][1(222t t t t t t q D tM t q M t K R q D t M t ∆-∆+∆-∆-∆--=∆+∆ (2-4)2.2中心差分法的特点(1)}{t t q ∆+的解是基于利用在时刻t 的平衡条件，即}{t t q ∆+是在假定t 时刻式(2-3)成立的条件下来计算的，该积分过程称为显式积分法，因此中心差分法为显式差分法。
但是，用式(2-4)计算}{t t q ∆+时，不仅要用到}{t q 而且还要用到}{t t q ∆-。
所以，计算在时刻t+△t 的解，必须存在一个具体的起始过程，即必须利用初始条件。
这样，已知}{0q 和}{0q，可由式(2-4)求得}{0q，由式(2-1)和式(2-3)可求得}{t q ∆-，如式(2-5)所示。
}{2}{}{}{0200q t q t q q t ∆+∆-=∆- (2-5) 我们称式(2-5)为差分格式(2-4)的初始条件。
(2)在求解时不需对刚度矩阵[K ]进行三角分解。
在应用此法时一般采用集束质量矩阵或称对角质量阵，而阻尼矩阵也通常为对角形式，这样，在运用式(2-4)时，就不需要对等号右端的系数矩阵t D t M ∆+∆2/][/]([2)进行三角分解，从而可以节省计算时间。
如果不考虑系统的阻尼，则0][=D ，这时式(2-4)可简化为：}{}]){[1(2t t t R q M t =∆∆- (2-6) 其中： }]{[1}]){[2]([}{}{22t t t t t q M t q M t K R R ∆-∆-∆--= (2-7) 从式(2-6)和式(2-7)可以看出，若质量阵[M ]为对角矩阵，即在系统中只考虑集束质量的情况，则方程(2-6)实际上是解耦的，即方程组中是相互独立的，因此只需进行矩阵相乘就可得到。
把式(2-6)改写为：}{][}{12t t t R M t q -∆-∆=写成分量形式，即：)(2)(i t ii i R m t q tt ∆=∆+ (2-8) 其中，)(i t t q ∆+和)(i t R 分别为}{t t q ∆+和}{t R 的第i 个分量，而ii m 为质量矩阵的第i 个对角线元素。
对于0=ii m ，我们可以用静凝聚的方法进行处理。
若对总体刚度和质量矩阵不必进行三角分角，就意味着不必形成总体的[K ]和[M ]，因此，式(2-7)的计算只需在单元一级进行即可，即：}){2}]({[1})]{([}{}{2t t t e e e e e t t q q M tq K R R -∆--=∆-∑∑ (2-9) 式中[K e ]和[M e ]分别为与待求节点自由度有关的单元刚度矩阵与质量矩阵。
所以，使用式(2-8)和式(2-9)构成的中心差分具有明显的优点。
它不需要计算总体刚度矩阵[K ]和质量矩阵[M ]，求解过程在单元一级进行，计算效率较高。
而且，如果所有单元刚度和质量矩阵都相同时，或有大部分的单元都相同时，可以用一个或少数单元矩阵来求解整个系统。
在工程中，常会遇到重复对称式结构，这时用这种方法来求解将具有很高的效率。
(3)加速度和速度的差分格式都具有(△t 2)阶的误差。
(4)可以证明，中心差分法所要求的步长为： πncr T t t ≤∆≤∆ (2-10)其中，T n 为有限无限无集合体的最小周期，n 为单元系统的阶，cr t ∆为临界时间步长。
周期T n 可利用如下任何一种方法来计算：(1) 广义Jacobi 法；(2) 逆迭代法；(3) 正迭代法；(4) Rayleigh 商迭代法；(5) 多项式迭代法；(6) 斯图姆序列对分法；(7) 行列式搜索法；(8) 子空间迭代法。
事实上，T n 为系统最大特征值的倒数。
要求使用的时间步长△t 小于临界时间步长cr t ∆的差分格式，称为条件稳定的。
若使用一个大于cr t ∆的时间步长，则积分将是不稳定的，此时由数值积分在计算机上导致的舍入误差会增大，从而使结果失去意义。
中心差分法是一种条件稳定的方法，这一缺点导致了它在使用上具有很大的局限性。
首选是时间步长选择困难，当矩阵阶数较高时，而且同时原系统的质量分布不均匀，将使得[M ]矩阵的个别对角上的元素很小，从而使系统的高阶特征趋于无限大，T n 趋于零，这样便使积分过程成为不可能。
所以，仅仅因为一个十分小的质量元素，就会导致计算效率的降低。
当这种情况出现时，通常采用静凝聚的方法使矩阵的阶数降低，但这仅仅能在低频段模拟原系统。
因此，中心差分法是条件稳定的非线性动力学分析方法，有其局限性。
3、Wilson -θ法Wilson -θ法是线性加速度方法的扩展，基于以下两条假定：(1)当时间从t 至变化到时间(t t ∆+θ)时，加速度从t q变化为t t q ∆+θ ； (2)在时间区段t ∆θ内，结构的刚度、阻尼、地面运动加速度都无改变。
则由上述两条假定可得：}){}({}{}{t t t t t q q tq q -∆+=∆++θτθτ (3-1) 式中，θ为参数，可根据积分的稳定性和精度进行选取，一般要求大于1。
对式(3-1)两端的τ进行积分，可得：}){}({2}{}{}{2t t t t t t q q tq q q -∆++=∆++θτθττ (3-2) 再对式(3-2)两端的τ进行积分，可得：}){}({6}{21}{}{}{22t t t t t t t q q tq q q q -∆+++=∆++θτθτττ (3-3) 在式(3-2)和式(3-3)中令t ∆=θτ，则可得在(t t ∆+θ)时刻的速度和位移分别如式(3-4)和式(3-5)所示：}){}({2}{}{t t t t t t q qtq q +∆+=∆+∆+θθθ (3-4) }){2}({6}{}{}{22t t t t t t t q qt t q q q +∆+∆+=∆+∆+θθθθ (3-5) 将加速度}{t t q∆+θ 和速度}{t t q ∆+θ 写成由位移}{t t q ∆+θ表达的形式，如式(3-6)和式(3-7)所示：}{2}{6}){}({6}{22t t t t t t t q q t q q t q -∆--∆=∆+∆+θθθθ (3-6) }{2}{2}{}({3}{t t t t t t t q t q q q t q ∆---∆=∆+∆+θθθθ (3-7) Wilson -θ法认为在(t t ∆+θ)时刻，系统满足动力平衡方程(1-1)，并由Wilson -θ法的第二条假定，我们把式(3-6)和式(3-7)代入到式(1-1)中，可得到式(3-8)。
[]}{}]{[}]{[}{t t t t t t t t R q K q D qM ∆+∆+∆+∆+=++θθθθ (3-8) 其中，}){}({}{}{t t t t t t R R R R -+=∆+∆+θθ (3-9)求(t t ∆+θ)时刻的位移、速度和加速度是我们所要得到的。
求解的过程如下所示：(1) 将式(3-6)和式(3-7)代入到式(3-8)中，并加以整理可得：}]){[3][)(6]([2t t q D tM t K ∆+∆+∆+θθθ =}){2}{6}{)(6]([}){}({}{2t t t t t t t q q t q t M R R R +∆+∆+-+∆+θθθ +}){2}{2}{3]([t t t q t q q t D ∆++∆θθ (3-10) 由式(3-10)可以求得}{t t q ∆+θ。
(2) 将式(3-6)分别代入式(3-1)、式(3-2)和式(3-3)中，并令t ∆=τ，可得：}){31(}{)(6}){}({)(6}{22tt t t t t t q q t q q t q θθθθ-+∆--∆=∆+∆+ (3-11) }){}({2}{}{t t t t t t q q t q q +∆+=∆+∆+ (3-12) }){2}({6}{}{}{2t t t t t t t q q t q t q q +∆+∆+=∆+∆+ (3-13) Wilson -θ法的特点：(1) 该法是一种隐式积分法，因为刚度矩阵[K ]是未知位移矢量的系数矩阵。