差分作业二

数学建模作业三——水位问题：

问题：

某居民小区有一个直径10m 的圆柱形水塔，每天午夜24时向水塔供水，此后每隔2h 记录水位，如下表，计算小区在这些时刻每小时的用水量。

模型及其求解：

令()f x 为x 时刻水塔的水量，根据题目问题，则可以先求出每个时刻的水量差即()f x —(1)f x +。

根据Taylor theorem ： ()''()()(1)()'()...()2!!

n n f x f x f x f x f x o x n +-=++++ 又有 (3)''()()(1)()'()2!3!

f x f f x f x f x ξ+-=++，224ξ<< 则先求出'()f x ，''()f x 即可得出的每个时刻每小时水量差，结果较为准确，其

误差范围可控制在(3)(3)min[()]max[()],3!3!f x f x ?? ???

中。

根据离散的数据求解'()f x ，非两端的数据利用中点公式，即：

()()'()2f a h f a h f a h

+--≈

而两端的数据则利用三点公式，即：

0122103443'(),'()22n n n n y y y y y y f x f x h h

---+--+≈

≈ 同理可求得''()f x 的值。

根据题目所给数据，利用MATLAB计算并作图，程序如下：

function y=cxy2(x0,f)

g=length(x0)-2;

m=length(x0);

x=0.25*pi*x0;

for k=1:g

r(k+1)=-(x(k+2)/(f*2)-x(k)/(f*2));

end

r(1)=-((-3*x(1)+4*x(2)-x(3))/(2*f));

r(12)=-((x(10)-4*x(11)+3*x(12))/(2*f));

y=r(1:m)';

k=(2:2:24)';

y=cxy2([305 298 290 265 246 225 207 189 165 148 130 114],2);

a=y;

N=[k,a]

f x的负值：

得到各个时刻'()

再根据将值代回函数cxy2中：

y=cxy2([3.25 3.75 8.25 11 10 9.75 9 10.5 10.25 8.75 8.5 7.5],2);

b=y;

M=[k,b]

f x的负值，数据如下：

得出各个时刻''()

再根据

''() (1)()'()

2!

f x

f x f x f x

+-≈+

又有

''()

(1)()'()

2!

f x

f x f x f x

+-≈--，可求出各个时刻每小时用水量的修正值。

利用MATLAB计算并作图，程序如下：

k=(2:2:24)';

a=[305 298 290 265 246 225 207 189 165 148 130 114];

y=cxy2([305 298 290 265 246 225 207 189 165 148 130 114],2);

b=y;

y=cxy2([3.25 3.75 8.25 11 10 9.75 9 10.5 10.25 8.75 8.5 7.5],2);

c=y;

d=b+c/2 ;

H=[k,d]

subplot(2,2,1),plot(k,a,'*',k,a),grid,xlabel('time/h'),ylabel('water

level/cm'),title('original data')

subplot(2,2,2),plot(k,b,'*',k,b),grid,xlabel('time/h'),ylabel('water level/cm'),title('1S rank derivative')

subplot(2,2,3),plot(k,c,'*',k,c),grid,xlabel('time/h'),ylabel('water level/cm'),title('2S rank derivative')

subplot(2,2,4),plot(k,d,'*',k,d),grid,xlabel('time/h'),ylabel('water

level/cm'),title('correct graph of data')

得到的数据，以及原始数据、一阶导数负值、二阶导数负值、修正后的图形如下：

结论分析：

2h用水2.84713m，4h用水2.45443m，6h用水5.76783m8h用水8.46763m，10h用水7.97673m，12h用水7.75583m，14h用水6.99503m，16h用水8.12403m，18h用水8.22213m，20h用水7.04403m，22h用水6.79863m，24h用水6.16053m。

其中修正前4h用水2.94523m大于修正前2h用水量，而修正后4h用水2.45443m则小于修正后2h用水量，按生活常识凌晨2时用水量大于凌晨4时，是更合理的。

其中使用的中心公式与三点公式误差与泰勒定理误差可能会相互削弱，也可能相互加强，所以在此可忽略前者带来的误差，而后者的误差可以以高阶导数来

减小，所以在此我认为建立一个求导函数cxy2是有必要的，如本题中，(3)()

f x的

负值我们可以通过函数cxy2求得：

y=cxy2([0.5890 -0.9817 -1.4235 -0.3436 0.2454 0.1963 -0.1473 -0.2454 0.3436

0.3436 0.2454 0.5400],2);

c=y’/6

可得误差范围

(3)(3)

min[()]max[()]

,

3!3!

f x f x

??

?

??

为()

0.0546,0.1397

-，可知结果较

为精确。

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

中心差分法的基本理论与程序设计 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的函数形式。 中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为

有限差分法

有限差分法 一、单变量函数： 用中心差分法（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;%时间

偏微分中心差分格式实验报告(含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，计 算结果如下：

时域有限差分法(姚伟)介绍

伊犁师范学院硕士研究生 ————期末考核 科目：电磁波有限时域差分方法 姓名：姚伟 学号：1076411203009 学院：电子与信息工程学院 专业：无线电物理

时域有限差分法 1 选题背景 在多种可用的数值方法中，时域有限差分法(FDTD)是一种新近发展起来的可选方法。1966年，K.S.Yee 首次提出电磁场数值计算的新方法—时域有限差分法(Finite Difference- Time Domain ，简称FDTD)。经历了二十年的发展FDTD 法才逐渐走向成熟。上世纪80年代后期以来FDTD 法进入了一个新的发展阶段，即由成熟转为被广泛接受和应用的阶段。FDTD 法是解决复杂问题的有效方法之一，是一种直接基于时域电磁场微分方程的数值算法，它直接在时域将Maxwell 旋度方程用二阶精度的中心差分近似，从而将时域微分方程的求解转换为差分方程的迭代求解。是电磁场和电磁波运动规律和运动过程的计算机模拟。原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机内存容量要求较低、计算速度较快、尤其适用于并行算法。现在FDTD 法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量学、瞬态电磁场研究等多个领域[1]。 2 原理分析 2.1 FDTD 的Yee 元胞 E,H 场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理 t t ??=??=??E D H ε t t ??-=??-=??H B E μ 图1 Yee 模型 如图1所示，Yee 单元有以下特点[2]： 1）E 与H 分量在空间交叉放置，相互垂直；每一坐标平面上的E 分量四周由H 分量环绕，H 分量的四周由E 分量环绕；场分量均与坐标轴方向一致。 2）每一个Yee 元胞有8个节点，12条棱边，6个面。棱边上电场分量近似相等，用棱边的中心节点表示，平面上的磁场分量近似相等，用面的中心节点表示。 3）每一场分量自身相距一个空间步长，E 和H 相距半个空间步长 4）每一场分量自身相距一个时间步长，E 和H 相距半个时间步长，电场取n 时刻的值，磁场取n+0.5时刻的值；即：电场n 时刻的值由n-1时刻的值得到，磁场n+0.5时刻的值由n-0.5时刻的值得到；电场n 时刻的旋度对应n+0.5时刻的磁场值，磁场n+0.5时刻的

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

中心差分法计算单自由度体系的动力反应 北京工业大学结构工程 组员:胡建华 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)

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

中心差分法计算单自由度体系动力反应 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

差分方法基础

第二讲 有限差分法基本原理 一般的流体控制方程都是非线性的偏微分方程。在绝大多数情况下，这些偏微分方程无法得到精确解；而CFD 就是通过采用各种计算方法得到这些偏微分方程的数值解，或称近似解。当然这些近似解应该满足一定的精度。目前，主要采用的CFD 方法是有限差分法和有限体积法。本讲主要介绍有限差分法，它也是下一讲中的有限体积法的基础[1]。 有限差分法求解流动控制方程的基本过程是：首先将求解区域划分为差分网格，用有限个网格点代替连续的求解域，将待求解的流动变量（如密度、速度等）存储在各网格点上，并将偏微分方程中的微分项用相应的差商代替，从而将偏微分方程转化为代数形式的差分方程，得到含有离散点上的有限个未知变量的差分方程组。求出该差分方程组的解，也就得到了网格点上流动变量的数值解。 2.1 差分和逼近误差 由于通常数字计算机只能执行算术运算和逻辑运算，因此就需要一种用算术运算来处理函数微分运算的数值方法。而有限差分法就是用离散网格点上的函数值来近似导数的一种方法。 设有x 的解析函数)(x f y =，从微分学知道函数y 对x 的导数为 x x f x x f x y dx dy x x ?-?+=??=→?→?)()(lim lim 00 （2-1） dy 、dx 分别是函数及自变量的微分，dx dy /是函数对自变量的导数，又称微商。相应地，上式中的x ?、y ?分别称为自变量及函数的差分，x y ??/为函数对自变量的差商。在导数的定义中x ?是以任意方式逼近于零的，因而x ?是可正可负的。在差分方法中，x ?总是取某一小的正数。这样一来，与微分对应的差分可以有三种形式： 向前差分 )()(x f x x f y -?+=? 向后差分 )()(x x f x f y ?--=? 中心差分 )2 1()21(x x f x x f y ?--?+=? 上面谈的是一阶导数，对应的称为一阶差分。对一阶差分再作一阶差分，就得到二阶差分，记为y 2?。以前向差分为例，有

对流扩散方程有限差分方法

对流扩散方程有限差分方法 求解对流扩散方程的差分格式有很多种，在本节中将介绍以下3种有限差分格式：中心差分格式、Samarskii 格式、Crank-Nicolson 型隐式差分格式。 3.1 中心差分格式 时间导数用向前差商、空间导数用中心差商来逼近，那么就得到了（1）式的中心差分格式]6[ 2 1 11 1122h u u u v h u u a u u n j n j n j n j n j n j n j -+-+++-=-+-τ （3） 若令 h a τ λ=，2h v τ μ=，则（3）式可改写为 )2()(2 111111n j n j n j n j n j n j n j u u u u u u u -+-+++-+--=μλ （4） 从上式我们看到，在新的时间层1+n 上只包含了一个未知量1 +n j u ，它可以由时间层n 上的值n j u 1-，n j u ，n j u 1+直接计算出来。因此，中心差分格式是求解对 流扩散方程的显示格式。 假定),(t x u 是定解问题的充分光滑的解，将1 +n j u ，n j u 1+，n j u 1-分别在),(n j t x 处 进行Taylor 展开： )(),(),(211 ττO t u t x u t x u u n j n j n j n j +??? ?????+==++ )(2),(),(3 22211 h O x u h x u h t x u t x u u n j n j n j n j n j +? ???????+????????+==++ )(2),(),(3 22211 h O x u h x u h t x u t x u u n j n j n j n j n j +????????+????????-==-- 代入(4)式，有 2 1 11 1122),(h u u u v h u u a u u t x T n j n j n j n j n j n j n j n j -+-+++---+-= τ )()()(2222 h O v x u v h O a x u a O t u n j n j n j ?-????????-?+????????++????????=τ )()()(222h O v a O x u v x u a t u n j n j n j ?-++????????-??? ?????+????????=τ

中心差分法

(1)中心差分法求解反应谱程序： clc load('GM.txt'); %输入原始数据 n=length(GM); t=0.01; Tn=0.02:0.02:4; m=length(Tn); MaxAcc=zeros(1,m); MaxVel=zeros(1,m); MaxDis=zeros(1,m); Damp=0.05; %确定阻尼比 mm=1; for T=0.02:0.02:4 %中心差分法求解结构响应时程 Dis=zeros(1,n+1); Vel=zeros(1,n); Acc=zeros(1,n); w=2*pi/T; c=2*w*Damp; k=w*w; Acc(1,1)=-1*GM(1,2); Vel(1,1)=0; Dis(1,2)=0; Dis(1,1)=(t^2)/2*Acc(1,1); for i=2:1:n Ke=1/(t^2)+c/(2*t); P=-1*GM(i-1,2)-(1/(t*t)-c/(2*t))*Dis(1,i-1)-(k-2/(t*t))*Dis(1,i); Dis(1,i+1)=P/Ke; end for ii=2:1:n-1 Vel(ii)=(Dis(ii+2)-Dis(ii))/(2*t); Acc(ii)=(Dis(ii+2)-2*Dis(ii+1)+Dis(ii))/(t*t)+GM(ii); end MaxAcc(1,mm)=max(abs(Acc)); %求解确定周期下结构最大响应MaxVel(1,mm)=max(abs(Vel)); MaxDis(1,mm)=max(abs(Dis)); mm=mm+1; end figure(1) %绘制反应谱 plot(Tn,MaxAcc(1,:)) title('加速度反应谱') xlabel('自振周期(s)') ylabel('最大绝对加速度Sa(g)')

有限差分方法概述

有限差分法（Finite Difference Method，简称FDM）是数值方法中最经典的方法，也是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。 构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。 下面我们从有限差分方法的基本思想、技术要点、应用步骤三个方面来深入了解一下有限差分方法。 1.基本思想 有限差分算法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。在采用数值计算方法求解偏微分方程时，再将每一处导数由有限差分近似公式替代，从而把求解偏微分方程的问题转换成求解代数方程的问题，即所谓的有限差分法。 2.技术要点 如何根据问题的特点将定解区域作网格剖分；如何把原微分