中心差分法的基本理论与程序设计

中心差分法的基本理论与程序设计
中心差分法的基本理论与程序设计

中心差分法的基本理论与程序设计

1程序设计的目的与意义

该程序通过用C语言(部分C++语言)编写了有限元中用于求解动力学问题的中心差分法,巩固和掌握了中心差分法的基本概念,提高了实际动手能力,并通过实际编程实现了中心差分法在求解某些动力学问题中的运用,加深了对该方法的理解和掌握。

2程序功能及特点

该程序采用C语言(部分C++语言)实现了用于求解动力学问题的中心差分法,可以求解得到运动方程的解答,包括位移,速度和加速度。计算简便且在算法稳定的条件下,精度较高。

3中心差分法的基本理论

在动力学问题中,系统的有限元求解方程(运动方程)如下所示:

阵和结点载荷向量,并分别由各自的单元矩阵和向量集成。

与静力学分析相比,在动力分析中,由于惯性力和阻尼力出现在平衡方程中,因此引入了质量矩阵和阻尼矩阵,最后得到的求解方程不是代数方程组,而是常微分方程组。常微分方程的求解方法可以分为两类,即直接积分法和振型叠加法。

中心差分法属于直接积分法,其对运动方程不进行方程形式的变换而直接进

行逐步数值积分。通常的直接积分是基于两个概念,

中心差分法的基本思路是用有限差分代替位移对时间的求导,将运动方程中的速度和加速度用位移的某种组合表示,然后将常微分方程组的求解问题转换为代数方程组的求解问题,并假设在每个小的时间间隔内满足运动方程,则可以求得每个时间间隔的递推公式,进而求得整个时程的反应。

在中心差分法中,加速度和速度可以用位移表示,即:

(2

t =

?(2t

=

?

而得到。为此将加速度和速度的表达式代入上式中,即可得到中心差分法的递推

公式:

个离散时间点的解的递推公式,这种数值积分方法又称为逐步积分法。需要指出

和速度的表达式可知:

2

a +

中心差分法避免了矩阵求逆的运算,是显式算法,且其为条件稳定算法,利

型载荷引起的波传播问题的求解,而对于结构动力学问题则不太合适。

4 中心差分法的有限元计算格式

利用中心差分法逐步求解运动方程的算法步骤如下所示: 1. 初始计算

(1)

(2)

(3)

ω=

(4)

(5)

(6)

2.

(1)

(2)

(3)

5 程序设计

5.1 程序流程

图1 程序流程图

各子程序主要功能为:

ArrayLU :LU Inverse :求矩阵的转置矩阵; ArrayMVector :矩阵和向量的乘法; LUSolve

5.2 输入数据及变量说明

5.2.1 输入数据

该程序的原始输入数据应包括三个部分: (1)

(2)

(3) 确定时间

步长,其中为了保证该算法的稳定性,需要满足

ω=

5.2.2 变量说明

该程序的各个变量含义如下: (1) num ,timeStep ,dt

num ——矩阵维度; timeStep ——时间步数; dt ——时间步长;

(2) M ,C ,K ,X ,V ,A ,P ,MM ,PT ,c0,c1,c2,c3

M ——质量矩阵; C ——阻尼矩阵; K ——刚度矩阵; X ——位移矩阵; V ——速度矩阵;

A ——加速度矩阵; P ——载荷向量; MM ——有效质量矩阵; PT —— c0,c1,c2,c3——积分常数;

6算例

6.1问题描述

应用本程序计算一个三自由度系统,它的运动方程是:

?

时,利用公式

,可以计算得到

6.2理论计算

6.2.1中心差分法(理论解)

(1)

则起步条件为

??

??

对于每一时间步长,需先计算有效载荷:

由上式得到的每一时间步长的位移结果如表1所示:

表1

0.00 0.00 0.00 0.03 0.13 0.36 0.79 1.46 2.37

3.42 0.00 0.03 0.19 0.58 1.26 2.24 3.43

4.69

5.84

6.77

0.40

1.48

2.97

4.52

5.82

6.71

7.22

7.51

7.85

8.45

(2)

按照相同的步骤,所得结果如下:

再计算下去,位移将继续增大,这是不稳定的典型表现。其原因是在条件稳定的可能得到有意义的结果。

6.2.2振型叠加法(精确解)

采用振型叠加法对上述问题进行计算,可以得到该问题的精确解。首先应求解的广义特征值问题为:

按照一般的线性代数方法可以得到上式的解答为:

利用上式,3个互不耦合的运动方程,即:

3

+

()

i t=

利用无阻尼情形的杜哈美积分公式,可以得到上述方程的解答为:

最后利用振型叠加得到系统的位移为:

根据上式计算得到的每一时间步长的位移值是系统响应的精确解,如表2 和表3所示。

(1)

表2

0.00 0.00 0.01 0.04 0.14 0.37 0.79 1.45 2.32 3.34 0.00 0.04 0.21 0.59 1.25 2.21 3.38 4.63 5.80 6.76 0.39

1.45

2.92

4.48

5.81

6.74

7.29

7.60

7.92

8.46

(2)

表3

3.89 3.46 -1.19 2.25 1.83 -2.40 1.78 2.25 -1.23 3.40 6.75 6.76 0.00 6.73 6.77 0.00 6.72 6.79 0.00 6.70 8.17

8.41

0.60

8.97

9.24

1.20

9.19

9.05

0.61

8.36

6.3 程序计算

(1)

中心差分法的基本理论与程序设计

中心差分法的基本理论与程序设计 1程序设计的目的与意义 该程序通过用C语言(部分C++语言)编写了有限元中用于求解动力学问题的中心差分法,巩固和掌握了中心差分法的基本概念,提高了实际动手能力,并通过实际编程实现了中心差分法在求解某些动力学问题中的运用,加深了对该方法的理解和掌握。 2程序功能及特点 该程序采用C语言(部分C++语言)实现了用于求解动力学问题的中心差分法,可以求解得到运动方程的解答,包括位移,速度和加速度。计算简便且在算法稳定的条件下,精度较高。 3中心差分法的基本理论 在动力学问题中,系统的有限元求解方程(运动方程)如下所示: ()()()() Ma t Ca t Ka t Q t ++= 式中,() a t分别是系统的结点加速度向 a t是系统结点位移向量,() a t和() 量和结点速度向量,,, M C K和() Q t分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量,并分别由各自的单元矩阵和向量集成。 与静力学分析相比,在动力分析中,由于惯性力和阻尼力出现在平衡方程中,因此引入了质量矩阵和阻尼矩阵,最后得到的求解方程不是代数方程组,而是常微分方程组。常微分方程的求解方法可以分为两类,即直接积分法和振型叠加法。 中心差分法属于直接积分法,其对运动方程不进行方程形式的变换而直接进行逐步数值积分。通常的直接积分是基于两个概念,一是将在求解域0t T内的任何时刻t都应满足运动方程的要求,代之仅在一定条件下近似地满足运动方程,例如可以仅在相隔t?的离散的时间点满足运动方程;二是在一定数目的t?区域内,假设位移a、速度a、加速度a的函数形式。 中心差分法的基本思路是用有限差分代替位移对时间的求导,将运动方程中的速度和加速度用位移的某种组合表示,然后将常微分方程组的求解问题转换为

有限差分法及其应用

有限差分法及其应用 1有限差分法简介 有限差分法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。该方程将解域划分为差分网格,用有限个网络节点代替连续的求解域。有限差分法通过泰勒级数展开等方法,把控制方程中的导数用网格节点上的函数值得差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。 2有限差分法的数学基础 有限差分法的数学基础是用差分代替微分,用差商代替微商而用差商代替微商的意义是用函数在某区域内的平均变化率来代替函数的真是变化率。而根据泰勒级数展开可以看出,用差商代替微商必然会带来阶段误差,相应的用差分方程代替微分方程也会带来误差,因此,在应用有限差分法进行计算的时候,必须注意差分方程的形式,建立方法及由此产生的误差。 3有限差分解题基本步骤 有限差分法的主要解题步骤如下: 1)建立微分方程 根据问题的性质选择计算区域,建立微分方程式,写出初始条件和边界条件。 2)构建差分格式 首先对求解域进行离散化,确定计算节点,选择网格布局,差分形式和步长;然后以有限差分代替无线微分,以差商代替微商,以差分方程代替微分方程及边界条件。 3)求解差分方程 差分方程通常是一组数量较多的线性代数方程,其求解方法主要包括两种:精确法和近似法。其中精确法又称直接发,主要包括矩阵法,高斯消元法及主元素消元法等;近似法又称间接法,以迭代法为主,主要包括直接迭代法,间接迭代法以及超松弛迭代法。4)精度分析和检验 对所得到的数值进行精度与收敛性分析和检验。 4商用有限差分软件简介 商用有限差分软件主要包括FLAC、UDEC/3DEC和PFC程序,其中,FLAC是一个基于显式有限差分法的连续介质程序,主要用来进行土质、岩石和其他材料的三维结构受力特性模拟和塑性流动分析;UDEC/3DEC是针对岩体不连续问题开发,用于模拟非连续介质在静,动态载荷作用下的反应;PFC是利用显式差分算法和离散元理论开发的微、细观力学程序,它是从介质的基本粒子结构的角度考虑介质的基本力学特性,并认为给定介质在不同应力条件下的基本特征主要取决于粒子之间接粗状态的变化,适用于研究粒状集合体的破裂和破裂发展问题,以及颗粒的流动(大位移)问题。

《有限差分法在微分方程中的应用》课程论文

课程论文

有限差分法在微分方程中的应用 本学期学习了《微分方程数值解》,本书中有限差分法给我留下的印象比较深刻,下边说说自己在方面的一点理解,请老师指正。 1.有限差分法的基本思想: 当系统的数学模型建立后,我们面对的主要问题就是微分积分方程的求解。基本思想是用离散的只含有限个未知量的差分方程组去近似地代替连续变量的微分方程和定解条件,并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似,对于方程中的积分用求和或及机械求积公式来近似代替,从而把原微分积分方程和边界条件转化成差分方程组。 2.有限差分法求解偏微分方程的步骤: 区域离散,即把所给偏微分方程的求解区域细分成由有限个格点组成的网格,这些离散点称作网格的节点; 近似替代,即采用有限差分公式替代每一个格点的导数。 逼近求解,换而言之,这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程。 从原则上说,这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格,或者通过离散点上的函数值进行插值计算来近似得到。理论上,当网格步长趋近于零时,差分方程组的解应该收敛于精确解,但由于机器字节的限制,网格步长不可能也没有必要取得无限小,那么差分法的收敛性或者说算法的稳定性就显得至关重要。因此,在运用有限差分法时,除了要保证精度外,还必须要保证其收敛性。 3.构造差分法的几种形式: 主要草用的是泰勒级数展开的方法。其基本差分表达式主要有三种形式:一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等。其中前两种形式为一阶计算精度,后一种为二阶计算精度。

笔记:线性常差分方程基本知识

本材料是关于线性常差分方程基本知识的笔记,参考了两个文献: 1、《差分方程》【日】福田武雄著穆鸿基译上海科学技术出版社1962年9月第一版 2、《常差分方程》王联、王慕秋著新疆大学出版社1991年2月第一版

目录 第一节差分 第二节和分 第三节对步长及定义域的约定 第四节阶乘多项式与差分 第五节Bernoulli多项式与差分 第六节几个公式,例题 第七节n阶线性常差分方程的解的结构 第八节 Lagrange变易常数法 第九节解n阶常系数齐次线性方程的特征根方法 第十节常系数对称型线性方程的解 第十一节几种特殊常系数非齐次线性方程的解法

第一节 差分 定义1.1:设函数()x f 的定义域是D ,R D ?,R x ∈?,0≠?x ,D x ∈?有D x x ∈?+,定义算子?为 ()()()x f x x f x f -?+=? 称x ?是x 的变化步长,()x f ?是()x f 在x 处的步长为x ?的一阶差分、阶差、有限差;D x ∈,函数()x f ?称为D 上的差分函数,简称差分;算子?是步长为x ?的差分算子。定义为 ()()x x f x f ?+=E 称()x f E 是()x f 在x 处的步长为x ?的一阶位移;称函数()x f E 是D 上的位移函数,简称位移;算子E 是步长为x ?的位移算子。定义算子I 为 ()()x f x f =I 称算子I 为恒等算子。称函数 ()x x f ??是D 上的差商函数,简称差商。 约定算子?与算子E 的步长相等。 注1.1: 大写希腊字母?、E 、I 的小写形式是δ、ε、ι,其英文单词形式是delta /`delt ?/ 、epsilon /ep`sail ?n/ 、 iota /ai`?ut ?/ 。 若D x ∈?,有D x x ∈?+,则N n ∈?,有D x n x ∈?+。 定理1.1:算子?、E 、I 有以下关系: ①()()()()()x f x f x f x f I -E =I -E =?,即I -E =?。 ②()()()()()x f x f x f x f I +?=I +?=E ,即I +?=E 。 ③()()()()x f x f E ?=?E ,即?E =E?。 定理1.2:算子?、E 是线性算子。对R b a ∈,,函数()x f 与()x g ,有以下等式 ()()()()()x g b x f a x bg x af ?+?=+? ()()()()()x g b x f a x bg x af E +E =+E 定义1.2:设N n ∈,作递推定义 ()()()x f x f x f =I =?0,()()() x f x f n n ??=?+1

有限差分法

有限差分法 一、单变量函数: 用中心差分法(matlab程序见附录)计算结果如下: 图1 中心差分法

表1 数据对比 二、一维热传导: 在此取φ(x)=0,g1(t)= g2(t)=100-100*exp(-t)

问题描述: 已知厚度为l的无限大平板,初温0度,初始瞬间将其放于温度为100度的流体中,流体与板面间的表面传热系数为一常数。 试确定在非稳态过程中板内的温度分布。 (1)显式差分法: 图3 显式差分法 (2)隐式差分法: 图4 隐式差分法

小结:显式格式仅当时格式是稳定的。(其中称为网格比) 隐式格式从k层求k+1层时,需要求解一个阶方程组。而且隐式格式的稳定性对网格比没有要求,即为绝对稳定的。 三、Possion方程: 取f=1,R=1 图5差分法

图6 误差小结:观察误差曲面,其绝对误差数量级为

附Matlab程序: 第1题: %===========================Boundary Value Problem 1 clear;clc; A=[-2.01 1 0 0 0 0 0 0 0; 1 -2.01 1 0 0 0 0 0 0; 0 1 -2.01 1 0 0 0 0 0; 0 0 1 -2.01 1 0 0 0 0; 0 0 0 1 -2.01 1 0 0 0; 0 0 0 0 1 -2.01 1 0 0; 0 0 0 0 0 1 -2.01 1 0; 0 0 0 0 0 0 1 -2.01 1; 0 0 0 0 0 0 0 1 -2.01;]; c1=[0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9]; C=0.01*c1-1*[0;0;0;0;0;0;0;0;1]; y=A\C; x=0:0.1:1; yn=[0;y;1]; ye=2*(exp(x)-exp(-x))/(exp(1)-exp(-1))-x; figure(1); plot(x,yn,'*',x,ye); legend('numerical solution','exact solution') xlabel('x','fontsize',20); ylabel('y','fontsize',20); set(gca,'fontsize',18); figure(2); err=abs(ye'-yn); plot(x,err); legend('error') xlabel('x','fontsize',20); ylabel('y','fontsize',20); set(gca,'fontsize',18); 第2题: %========================Boundary Value Problem 1_Explicit %显式 clear;clc l=20;%板厚 h=1;%步长 J=l/h; T=50;%时间

有限差分法

利用有限差分法分析电磁场边界问题 在一个电磁系统中,电场和磁场的计算对于完成该系统的有效设计师极端重要的。例如,在系统中,用一种绝缘材料是导体相互隔离是,就要保证电场强度低于绝缘介质的击穿强度。在磁力开关中,所要求的磁场强弱,应能产生足够大的力来驱动开关。在发射系统中进行天线的有效设计时,关于天线周围介质中电磁场分布的知识显然有实质性的意义。 为了分析电磁场,我们可以从问题所涉及的数学公式入手。依据电磁系统的特性,拉普拉斯方程和泊松方程只能适合于描述静态和准静态(低频)运行条件下的情况。但是,在高频应用中,则必须在时域或频域中求解波动方程,以做到准确地预测电场和磁场,在任何情况下,满足边界条件的一个或多个偏微分方程的解,因此,计算电池系统内部和周围的电场和磁场都是必要的。 对电磁场理论而言,计算电磁场可以为其研究提供进行复杂的数值及解析运算的方法,手段和计算结果;而电磁场理论则为计算电磁场问题提供了电磁规律,数学方程,进而验证计算结果。常用的计算电磁场边值问题的方法主要有两大类,其每一类又包含若干种方法,第一类是解析法;第二类是数值法。对于那些具有最简单的边界条件和几何形状规则的(如矩形、圆形等)问题,可用分离变量法和镜像法求电磁场边值问题的解析解(精确解),但是在许多实际问题中往往由于边界条件过于复杂而无法求得解析解。在这种情况下,一般借助于数值法求解电磁场的数值解。 有限差分法,微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网络来代替,这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 差分运算的基本概念: 有限差分法是指用差分来近似取代微分,从而将微分方程离散成为差分方程组。于是求解边值问题即转换成为求解矩阵方程[5]。 对单元函数 ()x f而言,取变量x的一个增量x?=h,则函数()x f的增量可以表示为 ()x f? = ()h x f+-()x f 称为函数()x f 的差分或一阶差分。函数增量还经常表示为 ()x f? = ? ? ? ? ? + 2 h x f - ? ? ? ? ? - 2 h x f

第二章计算流体力学的基本知识

第二章计算流体力学的基本知识 流体流动现象大量存在于自然界及多种工程领域中,所有这些工程都受质量守恒、动量守恒和能量守恒等基本物理定律的支配。这章将首先介绍流体动力学的发展和流体力学中几个重要守恒定律及其数学表达式,最后介绍几种常用的商业软件。 2.1计算流体力学简介 2.1.1计算流体力学的发展 流体力学的基本方程组非常复杂,在考虑粘性作用时更是如此,如果不靠计算机,就只能对比较简单的情形或简化后的欧拉方程或N-S方程进行计算。20 世纪30~40 年代,对于复杂而又特别重要的流体力学问题,曾组织过人力用几个月甚至几年的时间做数值计算,比如圆锥做超声速飞行时周围的无粘流场就从1943 年一直算到1947 年。 数学的发展,计算机的不断进步,以及流体力学各种计算方法的发明,使许多原来无法用理论分析求解的复杂流体力学问题有了求得数值解的可能性,这又促进了流体力学计算方法的发展,并形成了"计算流体力学" 。 从20 世纪60 年代起,在飞行器和其他涉及流体运动的课题中,经常采用电子计算机做数值模拟,这可以和物理实验相辅相成。数值模拟和实验模拟相互配合,使科学技术的研究和工程设计的速度加快,并节省开支。数值计算方法最近发展很快,其重要性与日俱增。 自然界存在着大量复杂的流动现象,随着人类认识的深入,人们开始利用流动规律来改造自然界。最典型的例子是人类利用空气对运动中的机翼产生升力的机理发明了飞机。航空技术的发展强烈推动了流体力学的迅速发展。 流体运动的规律由一组控制方程描述。计算机没有发明前,流体力学家们在对方程经过大量简化后能够得到一些线形问题解读解。但实际的流动问题大都是复杂的强非线形问题,无法求得精确的解读解。计算机的出现以及计算技术的迅速发展使人们直接求解控制方程组的梦想逐步得到实现,从而催生了计算流体力

微分方程与差分方程详细讲解与例题

第七章 常微分方程与差分方程 常微分方程是高等数学中理论性和应用性都较强的一部分,是描述客观规律的一种重要方法,是处理物理、力学、几何等应用问题的一个重要工具,微分和积分的知识是研究微分方程的基础。微分方程作为考试的重点容,每年研究生考试均会考到。特别是微分方程的应用问题,既是重点,也是难点,在复习时必须有所突破。 【数学一大纲容】常微分方程的基本概念;变量可分离的方程;齐次方程;一阶线性方程;伯努利(Bernoulli )方程;全微分方程;可用简单的变量代换求解的某些微分方程;可降阶的高阶微分方程;线性微分方程解的性质及解的结构定理;二阶常系数齐次线性微分方程;高于二阶的某些常系数齐次线性微分方程;简单的二阶常系数非齐次线性微分方程;欧拉(Euler )方程;微分方程的简单应用。 【数学二大纲容】常微分方程的基本概念;变量可分离的方程;齐次方程;一阶线性微分方程;可降阶的高阶微分方程;线性微分方程解的性质及解的结构定理;二阶常数齐次线性微分方程;高于二阶的某些常系数齐次线性微分方程;简单的二阶常系数非齐次线性微分方程;微分方程的一些简单应用。 【大纲要求】要理解微分方程的有关概念,如阶、解、通解、特解、定解条件等,掌握几类方程的解法:如变量可分离方程,齐次方程,一阶线性微分方程,伯努利方程,可降阶方程等。理解线性微分方程解的性质和解的结构,掌握求解常系数齐次线性方程的方法,掌握求解某些自由项的常系数非齐次线性方程的待定系数法。了解欧拉方程的概念,会求简单的欧拉方程。会用微分方程处理物理、力学、几何中的简单问题。 【考点分析】本章包括三个重点容: 1.常见的一阶、二阶微分方程求通解或特解。求解常微分方程重要的是判断方程为哪种类型,并记住解法的推导过程。 2.微分方程的应用问题,这是一个难点,也是重点。利用微分方程解决实际问题时,若是几何问题,要根据问题的几何特性建立微分方程。若是物理问题,要根据某些物理定律建立微分方程,也有些问题要利用微元法建立微分方程。 3.数学三要求掌握一阶常系数线性差分方程的求解方法,了解差分与差分方程及其通解与特解等概念,会用差分方程求解简单的经济应用问题。 【考点八十三】形如()()y f x g y '=的一阶微分方程称为变量可分离微分方程。可分离变量的微分方程的解题程序: 当()0,()()()() dy g y y f x g y f x dx g y '≠=? =时,然后左、右两端积分 (),()dy f x dx C g y =+?? 上式即为变量可分离微分方程的通解。其中,C 为任意常数,1 ()() dy g y g y ? 表示函数的一个原函数,()f x dx ?表示函数()f x 的一个原函数. 【例7.1】微分方程1+++='y x xy y 的通解为____________。

偏微分中心差分格式实验报告(含matlab程序)

二阶常微分方程的中心差分求解 学校:中国石油大学(华东)理学院 姓名:张道德 一、 实验目的 1、 构造二阶常微分边值问题: 22,(),(), d u Lu qu f a x b dx u a u b αβ?=-+=<

11122 222222333222122112 100121012010012 00N N N u f q h h u f q h h h u f q h h h q u f h h ---???? ??+-???? ??? ???? ???????-+-? ?????? ???????????=-+? ?????? ???????????-???? ????????-+????? ?? ????? 可以看出系数矩阵为三对角矩阵,而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解,则可以得出二阶常微分方程问题的数值解。 四、 举例求解 我们选取的二阶常微分方程边值问题为: 2 22242,01 (0)1,(1), x d u Lu x u e x dx u u e ?=-+=-<

差分法

第三章 有限差分法 函数()f x ,x 为定义在区间[]a b ,上的连续变 量。将区间[]a b ,等分成n 份,令()h b a n =-称为 步长,x 在这些离散点处的取值为 x a ih i =+ ()i n =01,,, 称为节点。函数()f x 在这些节点处的差值 ()()()()()() f x h f x f x f x h f x h f x h i i i i i i +---+--??? ?? (5-1) 分别称为一阶向前、向后和中心差分,可以用它 们作为函数()f x 在x i 处的微分近似值。这些差分 与相应x 区间的比值 ()()[] ()() [ ] ()()[] 1 1 1 2h f x h f x h f x f x h h f x h f x h i i i i i i +---+--?????? ??? (5-2) 分别称为一阶向前、向后和中心差商,可以用它 们作为函数()f x 在x i 处的导数近似值。完全类似 地可以定义高阶差商,例如常用的二阶中心差商 ()()()[] 1 22h f x h f x f x h i i i +-+- (5-3) 可以作为函数()f x 在x i 处的二阶导数近似值。 §3.1 常微分方程初值问题的差分解法 考虑电学中的一个问题:如图5-1。研究 电容器上的电荷随时间的变化规律。 图5-1 RC 放电回路 这个问题对应的微分方程及其定解条件为:

d d Q t Q RC Q Q t =-=??? ??=00 (5-4) 这是一阶微分方程的初值问题,它的解析解为 Q Q e t RC =-0 (5-5) 一、欧拉(Euler )折线法 求解下列普遍形式的一阶微分方程的初值 问题: ()[]()'=∈=?????y f x y x a b y a y ,,0 (5-6) 首先,将区间[]a b ,等分n 份,取值 a x x x b n =<<<=01 ,步长h x x i i =-+1。 然后,用一阶向前差商近似一阶导数,即 ()() ()()[] y x y x h y x f x y x i i i i i +-≈'=1, (5-7) 简记()y x y i i ≈,则式(5-7)可以写成差分格式: ()y y h f x y i i i i +=+?1, ()i n =-011,,, (5-8) 此即向前欧拉差分格式。这是一个递推计算格式, 从区间左端点即式(5-6)中的初始条件出发,按式 (5-8)依次可以算到区间右端点,得到的 y y y n 12,,, 就是原方程解()y x 的近似值。 应用式(5-8)计算RC 放电方程(5-4),按SI 单 位制,取Q 010=,RC =8,时间步长h =1,计 算结果如下:

分歧理论及其应用

现代电路理论 -------分歧理论及其应用

分歧理论及其应用 引言:近二、三十年来,分歧现象(bifurcation phenomena)及理论(bifurcation theory)在数学及自然科学上受到格外的重视及研究。随着科学技术的迅速发展,非线性问题大量出现于自然科学、工程技术乃至社会科学的许多领域,成为当前科学研究的热点。分歧现象是普遍存在的,是非线性系统的重要特点之一,它普遍地存在于数学、物理学、化学、经济学、社会学、生态学等各个领域,像数学中的解不唯一、物理学中的相变、工程中的静力与动力失稳、经济学中的马太效应、电子学中的周期振荡等等,都可以从分歧的角度去研究[1]。 1.分歧理论概述 分歧理论是近半个世纪以来逐步形成的有重要应用价值的数学分支,它反映的是流的拓扑结构随参数的变化而引起的质的变异,不论在数学理论上还是在现实应用中都具有极为重要的意义。近半个世纪以来,分歧理论的研究一直受到人们的广泛关注,也得到了很大的发展。国际电力界从20世纪80年代开始研究和应用分歧理论,在电压稳定、轴系扭振以及低频振荡的研究中均取得了新的突破。在上个世纪七十年代初,Crandall和Rabinowitz的两个基本分歧定理是由隐函数定理证明的,至今在数学,生物,工程上广为应用[2]。 分歧的含义是:对于含参数的系统,当参数发生变动并经过某些临界值时,系统的定性性态(即其拓扑结构,例如平衡状态、解的数目、周期运动的数目以及稳定性等)发生突然变化的现象。从数学角度而言,分歧理论主要是研究非线性代数方程(微分方程、积分方程、差分方程等)中参数对解的定性性质的影响,其中参数与解的稳定性、周期性、平衡位置等基本性质的关系是研究重点。 2. 分歧的定义 首先我们来看看一个经常可见到的现象。拿一根细长的金属棒。在棒的两头向内稍稍用力,此时棒不会弯曲。当力量够大时,则棒会弯起来。再继续加大压力,棒可能会弯了两弯。其变化如下图:

差分方程及其应用

差分方程及其应用 在经济与管理及其它实际问题中,许多数据都是以等间隔时间周期统计的。例如,银行中的定期存款是按所设定的时间等间隔计息,外贸出口额按月统计,国民收入按年统计,产品的产量按月统计等等。这些量是变量,通常称这类变量为离散型变量。描述离散型变量之间的关系的数学模型成为离散型模型。对取值是离散化的经济变量,差分方程是研究他们之间变化规律的有效方法。 本章介绍差分方程的基本概念、解的基本定理及其解法,与微分方程的基本概念、解的基本定理及其解法非常类似,可对照微分方程的知识学习本章内容。 §1 基本概念 线性差分方程解的基本定理 一、 基本概念 1、函数的差分 对离散型变量,差分是一个重要概念。下面给出差分的定义。 设自变量t 取离散的等间隔整数值:,,,, 210±±=t t y 是t 的函数,记作)(t f y t =。显然,t y 的取值是一个序列。当自变量由t 改变到1+t 时,相应的函值之差称为函数 )(t f y t =在t 的一阶差分,记作t y ?,即 )()1(1t f t f y y y t t t -+=-=+?。 由于函数)(t f y t =的函数值是一个序列,按一阶差分的定义,差分就是序列的相邻值之差。当函数)(t f y t =的一阶差分为正值时,表明序列是增加的,而且其值越大,表明序列增加得越快;当一阶差分为负值时,表明序列是减少的。 例如:设某公司经营一种商品,第t 月初的库存量是)(t R ,第t 月调进和销出这种商品的数量分别是)(t P 和)(t Q ,则下月月初,即第1+t 月月初的库存量)1(+t R 应是 )()()()1(t Q t P t R t R -+=+, 若将上式写作 )()()()1(t Q t P t R t R -=-+, 则等式两端就是相邻两月库存量的改变量。若记 ))()1()(t R t R t R -+=?, 并将理解为库存量)(t R 是时间t 的函数,则称上式为库存量函数)(t R 在t 时刻(此处t 以月为单位)的差分。 按一阶差分的定义方式,我们可以定义函数的高阶差分。函数)(t f y t =在t 的一阶差

差分方程的基本知识(3)

差分方程模型的理论和方法 1、差分方程:差分方程反映的是关于离散变量的取值与变化规律。通过建立一个或几个离散变量取值所满足的平衡关系,从而建立差分方程。 差分方程就是针对要解决的目标,引入系统或过程中的离散变量,根据实际背景的规律、性质、平衡关系,建立离散变量所满足的平衡关系等式,从而建立差分方程。通过求出和分析方程的解,或者分析得到方程解的特别性质(平衡性、稳定性、渐近性、振动性、周期性等),从而把握这个离散变量的变化过程的规律,进一步再结合其他分析,得到原问题的解。 2、应用:差分方程模型有着广泛的应用。实际上,连续变量可以用离散变量来近似和逼近,从而微分方程模型就可以近似于某个差分方程模型。差分方程模型有着非常广泛的实际背景。在经济金融保险领域、生物种群的数量结构规律分析、疾病和病虫害的控制与防治、遗传规律的研究等许许多多的方面都有着非常重要的作用。可以这样讲,只要牵涉到关于变量的规律、性质,就可以适当地用差分方程模型来表现与分析求解。 3、差分方程建模:在实际建立差分方程模型时,往往要将变化过程进行划分,划分成若干时段,根据要解决问题的目标,对每个时段引入相应的变量或向量,然后通过适当假设,根据事物系统的实际变化规律和数量相互关系,建立每两个相邻时段或几个相邻时段或者相隔某几个时段的量之间的变化规律和运算关系(即用相应设定的变量进行四则运算或基本初等函数运算或取最运算等)等式(可以多个并且应当充分全面反映所有可能的关系),从而建立起差分方程。或者对事物系统进行划分,划分成若干子系统,在每个子系统中引入恰当的变量或向量,然后分析建立起子过程间的这种量的关系等式,从而建立起差分方程。在这里,过程时段或子系统的划分方式是非常非常重要的,应当结合已有的信息和分析条件,从多种可选方式中挑选易于分析、针对性强的划分,同时,对划分后的时段或子过程,引入哪些变量或向量都是至关重要的,要仔细分析、选择,尽量扩大对过程或系统的数量感知范围,包括对已有的、已知的若干量进行结合运算、取最运算等处理方式,目的是建立起简洁、深刻、易于求解分析的差分方程。在后面我们所举的实际例子中,这方面的内容应当重点体会。

编制中心差分法程序报告(结构工程研究生作业)

中心差分法计算单自由度体系的动力反应 北京工业大学结构工程 组员:胡建华 S201204111 马 恒 S201204112 陈相家S201204083 张力嘉S201204022 0前言 时域逐步积分法是数值分析方法,它只假设结构本构关系在一个微小的时间步距内是线性的。时域逐步积分法是结构动力问题中一个得到广泛研究的课题,并在结构动力反应计算中得到广泛应用。由于引进的假设条件不同,可以有各种不同的方法,比如中心差分法,线性加速度法,平均常加速度法,Wilson-θ法等,其中中心差分法精度最高。在本文中,通过编制中心差分法计算单自由度体系的动力反应的程序来了解其应用及稳定性。 1中心差分法原理 中心差分法的基本思路:是将运动方程中的速度向量和加速度向量用位移的某种组合来表示,将微分方程组的求解问题转化为代数方程组的求解问题,并在时间区间内求得每个微小时间区间的递推公式,进而求得整个时程的反应。 中心差分法只在相隔t ?一些离散的时间区间内满足运动方程,其基于有限差分代替位移对时间的求导(即速度和加速度),如果采用等时间步长,t t i ?=?,则速度与加速度的中心差分近似为: t u u u i i ?-= -+? 21 1 (a) 2 1 12t u u u u i i i ?+-= -+? ? (b) 而离散时间点的运动为 )(),(),(i i i i i i t u u t u u t u u ? ?? ?? ? === ( =i 0,1,2,3,……) 由体系运动方程为:0)()()(=++? ? ?t ku t u c t u m i (c) 将速度和加速度的差分近似公式(a )和式(b )代入式(c )可以得到i t 时刻的运动方程: 0221 12 11=+?-+?+--+-+i i i i i i ku t u u c t u u u m (d )

中心差分法计算程序编程.doc

中心差分法计算程序编程 姓名:张泽伟 学号: 电话: 一、中心差分法程序原理说明 1.1 中心差分法思路 中心差分法的基本思路:是将运动方程中的速度向量和加速度向量用位移的某种组合来表示,将微分方程组的求解问题转化为代数方程组的求解问题,并在时间区间内求得每个微小时间区间的递推公式,进而求得整个时程的反应。 1.2 中心差分法原理 中心差分法只在相隔t ?一些离散的时间区间内满足运动方程,其基于有限 差分代替位移对时间的求导(即速度和加速度),如果采用等时间步长,t t i ?=?, 则速度与加速度的中心差分近似为: t u u u i i ?-= -+?211 (a) 21 12t u u u u i i i ?+-= -+?? (b) 而离散时间点的运动为 ) (),(),(i i i i i i t u u t u u t u u ??????=== ( =i 0,1,2,3,……) 由体系运动方程为:0 )()()(=++???t ku t u c t u m i (c) 将速度和加速度的差分近似公式(a )和式(b )代入式(c )可以得到i t 时 刻的运动方程: 02211211=+?-+?+--+-+i i i i i i ku t u u c t u u u m (d ) 在(d )式中,假设i u 和1-i u 是已知的,即在i t 及i t 以前时刻的运动已知,则 可以把已知项移到方程的右边,整理得到: 12212)2()2()2(-+?-?-?--=?+?i i i u t c t m u t m k u t c t m (e)

由式(e )就可以根据i t 及i t 以前时刻的运动,求得1+i t 时刻的运动,如果需 要可以用式(a )和式(b )求得体系的速度和加速度。 1.3 初始条件转化 假设给定的初始条件为 ),0(), 0(00? ?==u u u u (g ) 由式(g )确定1-u 。在零时刻速度和加速度的中心差分公式为: t u u u ?-= -?2110 (h ) ` 21 0102t u u u u ?+-=-?? (i ) 将式(i )消去1u 得:020012???-?+?-=u t u t u u (j ) 而零时刻的加速度值0??u 可以用t =0时的运动方程 0000=++???ku u c u m 确定 即 )(1000ku u c m u --=?? ? (k ) 这样就可以根据初始条件 00,?u u 和初始荷载0P ,就可以根据上式确定1-u 的 值。 1.4 中心差分法编程思路 ① 基本数据准备和初始条件计算: )(1000ku u c m u --=?? ? 020012??? -?+?-=u t u t u u ② 计算等效刚度和中心差分计算公式中的相关系数: t c t m k ?+?=22

中心差分法在单自由度中的应用

中心差分法求解单自由度体系的自由振动问题 前言 时域逐步积分法是根据运动方程,引进某些假设,建立由t 时刻状态向量i u 、i u ?、i u ? ?到t +t ?时刻的状态向量1+i u 、1+?i u 、1+??i u 的递推关系,从而从t =0时刻的初始状态向量0u 、0?u 、0? ?u 出发,逐步求出各时刻的状态向量,由于引进的假设条件不同,可以有各种不同的方法,下面主要介绍一种时域逐步积分方法-中心差分法。 中心差分法(central difference method)原理[1] 中心差分法的基本思路:是将运动方程中的速度向量和加速度向量用位移的某种组合来表示,将微分方程组的求解问题转化为代数方程组的求解问题,并在时间区间内求得每个微小时间区间的递推公式,进而求得整个时程的反应。 中心差分法只在相隔t ?一些离散的时间区间内满足运动方程,其基于有限差分代替位移对时间的求导(即速度和加速度),如果采用等时间步长,t t i ?=?,则速度与加速度的中心差分近似为: t u u u i i ?-=-+?211 (a) 2 112t u u u u i i i ?+-=-+?? (b) 而离散时间点的运动为 )(),(),(i i i i i i t u u t u u t u u ??????=== ( =i 0,1,2,3,……) 由体系运动方程为:0)()()(=++???t ku t u c t u m i (c) 将速度和加速度的差分近似公式(a )和式(b )代入式(c )可以得到i t 时刻的运动方程: 02211211=+?-+?+--+-+i i i i i i ku t u u c t u u u m (d ) 在(d )式中,假设i u 和1-i u 是已知的,即在i t 及i t 以前时刻的运动已知,则可以把已知项移到方程的右边,整理得到: 12212)2()2()2(-+?-?-?--=?+?i i i u t c t m u t m k u t c t m (e)

有限差分法

有限差分法 有限差分法有限差分法 finite difference method 微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散 点构成的网格来代替,这些离散点称作网格的节点;把连续定解区域上的连续变量的函 数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差 分方程组,解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便 可以从离散解得到定解问题在整个区域上的近似解。 有限差分法的主要内容包括:如何根据问题的特点将定解区域作网格剖分;如何把原 微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和 计算结果的正确,还需从理论上分析差分方程组的性态,包括解的唯一性、存在性和差分 格式的相容性、收敛性和稳定性。对于一个微分方程建立的各种差分格式,为了有实用意义,一个基本要求是它们能够任意逼近微分方程,这就是相容性要求。另外,一个差分格 式是否有用,最终要看差分方程的精确解能否任意逼近微分方程的解,这就是收敛性的概念。此外,还有一个重要的概念必须考虑,即差分格式的稳定性。因为差分格式的计算过 程是逐层推进的,在计算第n+1层的近似值时要用到第n层的近似值,直到与初始值有关。前面各层若有舍入误差,必然影响到后面各层的值,如果误差的影响越来越大,以致 差分格式的精确解的面貌完全被掩盖,这种格式是不稳定的,相反如果误差的传播是可以 控制的,就认为格式是稳定的。只有在这种情形,差分格式在实际计算中的近似解才可能 任意逼近差分方程的精确解。关于差分格式的构造一般有以下3种方法。最常用的方法是 数值微分法,比如用差商代替微商等。另一方法叫积分插值法,因为在实际问题中得出的 微分方程常常反映物理上的某种守恒原理,一般可以通过积分形式来表示。此外还可以用 待定系数法构造一些精度较高的差分格式。 有限差分方法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。该方法 将求解域划分为差分网格,用有限个网格节点代替连续的求解域。有限差分法以Taylor 级数展开等方法,把控制方程中的导数用网格节点上的函数值的差商代替进行离散,从 而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数 问题的近似数值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。 对于有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式。从差分 的空间形式来考虑,可分为中心格式和逆风格式。 考虑时间因子的影响,差分格式还可以分为显格式、隐格式、显隐交替格式等。目 前常见的差分格式,主要是上述几种形式的组合,不同的组合构成不同的差分格式。差分 方法主要适用于有结构网格,网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。

差分方程在经济学中的应用应用数学

本科毕业论文(设计) 论文题目:差分方程在经济学中的应用 学生姓名:雷晶 学号: 1004970226 专业:数学与应用数学 班级:数学1002班 指导老师:舒蕊艳 完成日期:2014年5月20日

差分方程在经济学中的应用 内容摘要 本文叙述了研究差分方程的意义和背景、差分方程的定义、常见的解法以及差分方程相关模型,重点介绍差分方程经济学中的应用模型—筹措教育经费模型,包括问题的提出、模型举例和分析、提出假设、模型建立、模型求解、结果分析等等步骤对模型进行了更深层次的分析,做了进一步的推广. 本文所介绍的筹措教育经费模型主要研究的是子女的教育费用,假定某家庭从孩子m岁起,每月拿出一部分钱存进银行,用于投资子女的大学教育,并计划n年后支出一些,直到孩子大学毕业,全部用完账户中的资金. 差分方程的理论研究近十年来发展十分迅速,尤其是在经济领域,帮助人们解决了很多实际问题,筹措教育经费模型的建立为广大中国家庭子女教育的费用问题提供了明确的解决方法,是差分方程理论最贴近实际的模型之一. 关键词:差分方程存款模型经济增长模型筹措教育经费模型

, . . , , , , . a . ’s . , ’s ’s m n , . , . a . a ’s . 目录 一、绪论 (1) (一)研究差分方程在经济学中的应用的目的意义 (1) (二)研究背景 (2) 二、研究的理论基础 (2) (一)差分 (2) (二)差分方程 (3) (三)差分方程的解 (4) (四)特征根法 (4)

三、差分方程的经济应用模型简介 (5) (一)贷款模型 (5) (二)存款模型 (6) (三)乘数-加速数模型 (7) (四)哈罗德-多马经济增长模型 (10) (五)投入产出模型 (11) (六)筹措教育经费模型 (12) 四、总结 (14) 参考文献 (16)

中心差分法计算单自由度体系动力反应样本

中心差分法计算单自由度体系动力反应 1, 程序说明 中心差分法基于有限元查分代替位移对时间的求导( 即速度和加速度) 。如果采用等时间步长, i t t ?=?, 则速度和加速度的中心查分近似为 .11..11222i i i i i u u u t u u u u t +-+--=?-+=? 体系的运动方程为 ... ()()()()mu t cu t ku t P t ++= 联立以上三式, 得 1111222i i i i i i i u u u u u m c ku P t t +-+--+-++=?? 上式中, 假设i u 和1i u -是已知的, 即i t 和i t 以前时刻的运动已知, 则能够把已知 项移到方程的右边, 整理得 11222222i i i i m c m m c u P k u u t t t t t +-??????+=---- ? ? ???????????? 这样, 就能够计算体系任意时刻的位移, 速度和加速度。 2, 程序框图

3, 程序清单 %计算等效刚度和中心差分法计算公式中的系数 clear, m=17.5e3;k=875500;c=35000;aa=input(' 请选择时间步长 1 or 2 or 3 \n 1: dt=0.02 ; 2 : dt=0.3 ; 3: dt=其它\n'); if aa==1 dt=0.02; end if aa==2 dt=0.3; end if aa==3 dt=input('请输入时间步长\n dt= ') end t=0:dt:1.2; n=fix(1.2/dt+1);kr=m/(dt * dt) + c/(2 * dt); a=k-2 * m/(dt * dt); b=m/(dt * dt)-c/(2 * dt); %求力p p1=0:1:40;p2=39:-1:0;one=ones(1,40);p3=(one<0);p=1000*[p1,p2,p3]; for i=1:n if t(i)<=0.4,p(i)=100000*t(i); end if t(i)>0.4&&t(i)<=0.8,p(i)=8*t(i); end if t(i)>0.8,p(i)=0; end

相关文档
最新文档