【文档说明】计算机仿真教案02--第二章数值积分法的系统仿真-课件.ppt，共(65)页，946.519 KB，由小橙橙上传

转载请保留链接：https://www.ichengzhen.cn/view-76786.html

以下为本文档部分文字说明：

船舶与海洋工程学院第二章数值积分法的系统仿真仿真技术基础·数值积分法的系统仿真2.1概述2.1概述连续系统仿真，从本质上：对原连续系统从时间、数值两个方面对原系统进行离散化并选择合适的数值计算方法来近似积分运算在数值积分法的计算中，只计算了采样点的值，

相当于是对系统模型进行了离散化处理，所以从本质说，数值积分法也是离散化方法，只不过它是从数值积分的角度出发，没有明确提出“离散”这个概念仿真技术基础·数值积分法的系统仿真2.1概述1.相似原理其中u(t)为输入

变量，y(t)为系统变量；令仿真时时间隔为h，离散化后的输入变量为)(ˆntu系统变量为)(ˆnty其中nt表示t=nh),,(tuyfy设系统模型为如果)()(ˆnntutu)()(ˆnntyty即0)()(ˆ)(nnnu

tutute0)()(ˆ)(nnnytytyte（对所有n=0,1,2,…）则可认为两模型等价仿真技术基础·数值积分法的系统仿真2.1概述仿真技术基础·数值积分法的系统仿真2.1概述2.对仿真建模方法三个基本要求：（1）稳定性：不改变原系统的稳定性若原

连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的若原连续系统是不稳定的，则离散化后得到的仿真模型也应是不稳定的（2）准确性：有不同的准确性评价准则，最基本的准则是：绝对误差准则：)()(ˆ)(nnnyty

tyte相对误差准则：)(ˆ)()(ˆ)(nnnnytytytyte其中规定精度的误差量仿真技术基础·数值积分法的系统仿真2.1概述（3）快速性：若第n步计算所对应的系统时间间隔为3.数值积分算法：,1nnnt

th计算机由)(nty计算)(1nty需要的时间为Tn，若Tn=hn称为实时仿真Tnhn称为超实时仿真Tnhn称为亚实时仿真对),,(tuyfy，已知系统变量y的初始条件00)(yty求随时间变化的过程)(tyy－－初值问题仿真技术基础·数值积分法的系统仿真2.1概述计算过程：

由初始点00)(yty的)(00ytf，ttdtytfyty0),()(0f(t,y)f(t0,yo)ttt0t1仿真技术基础·数值积分法的系统仿真2.2数值积分法2.2数值积分法2.2.1欧拉法)],(,[)(tytfty00)0()(yytyt仿真技术基础·数值积分法的

系统仿真节点间距为步长，通常可采用等距节点，即取hi=h(常数)。2.2数值积分法要计算出解函数y(x)在一系列节点a=x0<x1<…<xn=b处的近似值)1,...,0(1nixxhiii1.欧拉公式：向前差商近似导数hxyxyxy)()()(01

0),()()()(000001yxfhyxyhxyxy1y记为)1,...,0(),(1niyxfhyyiiii计算yn+1时,只用到前一步的结果yn,因此属于单步法x0x1仿真技术基础·数值积分法

的系统仿真其实质如图所示,用折线来近似实际的曲线2.2数值积分法Euler方法的几何体现)(jjxyhy),(1jjjjyxhfyy仿真技术基础·数值积分法的系统仿真定义：若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度2.2数值积分法定义：在假设yi=y(xi)，

即第i步计算是精确的前提下，考虑的截断误差Ri=y(xi+1)yi+1称为局部截断误差/*localtruncationerror*/。欧拉法的局部截断误差：定义：欧拉法具有1阶精度211()[()()()][(,)

]iiiiiiiiRyxyyxhyxOhyhfxy2()Oh(),jjhyj设R为计算的求解公式第步的截断误差且1()()kkjjEhRh步的累计截断误差为该求解公式第则称khEk)(点上的总体截断误差即该求解公式在kx仿真技术基础·数值积分法的系统

仿真1.1引言-2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-0.2000.511.522.53210xxxx)()(1122xyyxyy1()Rh2()Rh3()Rh)(1hE)(2

hE)(3hE仿真技术基础·数值积分法的系统仿真（1）隐式欧拉法/*implicitEulermethod*/2.2数值积分法2.欧拉公式的改进：用向后差商近似导数由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式，

而前者称为显式/*explicit*/欧拉公式一般先用显式计算一个初值，再迭代求解隐式欧拉法的局部截断误差：即隐式欧拉公式具有1阶精度hxyxyxy)()()(011))(,()(1101xyxfhyxy)1,..

.,0(),(111niyxfhyyiiii211()()iiiRyxyOhx0x1仿真技术基础·数值积分法的系统仿真（3）中点欧拉公式/*midpointformula*/2.2数值积分法（2）梯形公式/*trapezoidformula—显、隐式两种算法的平均中心差商

近似导数)1,...,0()],(),([2111niyxfyxfhyyiiiiii注：即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。)()(311hOyxyRiiihxyxyx

y2)()()(021))(,(2)()(1102xyxfhxyxy1,...,1),(211niyxfhyyiiii假设，则可以导出即中点公式具有2阶精度。)(),(11iiiixyyxyy)()(311hOyxyRiii仿真

技术基础·数值积分法的系统仿真2.2数值积分法比较：方法显式欧拉隐式欧拉梯形公式中点公式简单精度低稳定性最好精度低,计算量大精度提高计算量大精度提高,显式多一个初值,可能影响精度仿真技术基础·数值积分法的系统仿真2.

2数值积分法3改进的欧拉法/*modifiedEuler’smethod*/Step2:再将代入隐式梯形公式的右边作校正，得到1iy)],(),([2111iiiiiiyxifyixifhyiy)1,...,0(),(,),(211niyxfhy

xfyxfhyyiiiiiiii此法亦称为预测-校正法/*predictor-correctormethod*/。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。它的稳定性高于显式欧拉法。Step1:

先用显式欧拉公式作预测，算出),(1iiiiyixifhyiy仿真技术基础·数值积分法的系统仿真收敛性/*Convergency*/2.2数值积分法4收敛性/*Convergency*/定义：若某算法对于任意固定的x=xi=x0+i

h，当h0(同时i)时有yiy(xi)，则称该算法是收敛的。例：就初值问题考察欧拉显式格式的收敛性。0)0(yyyy解：该问题的精确解为xeyxy0)(欧拉公式为iiiiyhyhyy)1(1对任意固定的x=xi=ih，有iixhhxihyhyy

])1[()1(/10/0ehhh/10)1(lim)(0ixxyeyi仿真技术基础·数值积分法的系统仿真2.2数值积分法例1.Euler1)0(102yxyxyy解:2(,)xfxyyy1,1,10,000y

bnax由前进Euler公式),(1jjjjyxhfyynj,,2,1)2(jjjjyxyhy0.1h用公式求初值问题显然取仿真技术基础·数值积分法的系统仿真2.2数值积分法得)2(00001yxyhyy1.1)1021(1.01)2(11112yx

yhyy1918.1)1.11.021.1(1.01.1],[yx01.00000.10001.10000.20001.19180.30001.27740.40001.35820.50001.43510.60001.50900.

70001.58030.80001.64980.90001.71781.00001.7848仿真技术基础·数值积分法的系统仿真2.2数值积分法用Euler公式的预测——校正系统求解例1例2解:由预测－

校正公式,有)2(00001yxyhyy1.1)1021(1.010110010122[()()]2hxxyyyyyy1.09181221121222[()()]2hxxyyyyyy)2(11112yxyhyy1827.1)0918.11.020918

.1(1.00918.11.1763仿真技术基础·数值积分法的系统仿真2.2数值积分法依此类推,得01.00000.10001.09180.20001.17630.30001.25460.40001.32780.5

0001.39640.60001.46090.70001.52160.80001.57860.90001.63211.00001.6819],[yx仿真技术基础·数值积分法的系统仿真考虑改进Euler法2.2数值积分法

2.2.2Runge-Kutta法),(111kkkkyxhfyy)],(),([2111kkkkkkyxfyxfhyy)(2211KKhyykk),(111kkyxfK),(112hKyxfKkk)(00xyy---------

-(1)仿真技术基础·数值积分法的系统仿真梯形公式具有2阶精度2.2数值积分法改进Euler法是由梯形公式和Euler公式复合而成同样可以证明,改进Euler法也具有2阶精度形如(1)式的求解公式称为二阶Runge-Kutta法对于Simpson求解公式：)],(

))2(,2(4),([611111kkkkkkkkyxfhxyhxfyxfhyy这是隐式多步法选取适当的显化方法,可得类似(1)的高阶Runge-Kutta方法以下使用中值定理进行

推导仿真技术基础·数值积分法的系统仿真对于常微分方程的初值问题2.2数值积分法一、Runge-Kutta方法的导出0)(),(yaybxayxfy的解),(xyy有上使用微分中值定理在区间,],[1nnxx))(()()(111nnnnnxxyxyx

y),(11nnnxx)()()(11nnnyhxyxy即仿真技术基础·数值积分法的系统仿真2.2数值积分法上的平均斜率在区间可以认为是],[)(1nnxxxyyKhKxyxynn)()(1----------(3)引入记号)

(1nyK)](,[11nnyfKxxxynn上平均斜率的近似值间在区出只要使用适当的方法求],[)(1就可得到相应的Runge-Kutta方法hKyynn1----------(4)即(4)式1

nxnxxy)(xyyK仿真技术基础·数值积分法的系统仿真2.2数值积分法二、低阶Runge-Kutta方法1nxnxxy)(xyy如下图Kxxxyxxynnn上的平均斜率在处的斜率作为在如果以],[)()(.111即)(1nxyK)]

(,[11nnxyxf则(4)式化为),(111nnnnyxhfyy),(11nnyxf即Euler方法Euler方法也称为一阶Runge-Kutta方法2()()nRhOh由于----(5)KK仿真技术基础·数值积分法的系统仿真2.2数值积分法1

nxnxxy)(xyy上的平均斜率在的算术平均值作为和处的斜率和在如果以],[)()(.21211nnnnxxxyKKxxxy)(11nxyK),(11nnyxf)(2nxyK)](,[nnxyxf),(nnyxf),(11hKyxfnn(由(

5)式)令221KKK则(4)式化为K1K2K仿真技术基础·数值积分法的系统仿真2.2数值积分法)(2211KKhyynn),(111nnyxfK),(112hKyxfKnn)(00xyy-----------(6)和(

1)式一致,即改进Euler公式,也称为二阶Runge-Kutta法3()()nRhOh仿真技术基础·数值积分法的系统仿真2.2数值积分法三、高阶Runge-Kutta方法2],[1212111hxxxxxnnnnn

上增加一点如果上的平均斜率在作为的加权平均值和、处的斜率和、在且以],[)()(1321211nnnnnxxxyKKKxxxxy)(11nxyK),(11nnyxf)(212nxyK)](,[2121nnxyxf)(3nxyK)](,[nnxyxf未

知K1K2K3K1nxnx21hxnxy)(xyy仿真技术基础·数值积分法的系统仿真2.2数值积分法)(21nxy)2,2(111Khyhxfnn令112Khyn)(212nxyK)2(121KKhyn)(nxy))2(,(1211KKhyhxfnn

)(3nxyK)(2111nnxyKx预测处的斜率如果以)(22111nnnxyKKhxx预测、处的斜率、同样以令参照Simpson求解公式Euler公式)],())2(,2(4),([61111

1nnnnnnnnyxfhxyhxfyxfhyy仿真技术基础·数值积分法的系统仿真2.2数值积分法取321616461KKKK则（4）式化为)4(63211KKKhyynn),(111nnyxfK)2,2(1112KhyhxfKnn)

)2(,(12113KKhyhxfKnn-----------(7))(00xyy(7)式称为三阶Runge-Kutta方法仿真技术基础·数值积分法的系统仿真2.2数值积分法展开式处的在、、分别作TailorxKKKn1321),(111nnyx

fK)(1nxy)2,2(1112KhyhxfKnn),(2),(2),(1111111nnynnxnnyxfhKyxfhyxf)()],(4),(4),(4[!213112121112112hOyxfKhy

xfKhyxfhnnynnxynnx)(2)(11nnxyhxy)()(8312hOxyhn仿真技术基础·数值积分法的系统仿真2.2数值积分法))2(,(12113KKhyhxf

Knn),()2(),(),(11121111nnynnxnnyxfKKhyxfhyxf)()],()2(),()2(),([!21311212211122112hOyxfKKhyxfKKhyxfhnnynnxynnx

),()22()()(111211nnynnyxfKKhxyhxy)()],()44(),()22()([2311212111212hOyxfKKKyx

fKKxyhnnynnxyn)()(2)()(31211hOxyhxyhxynnn仿真技术基础·数值积分法的系统仿真2.2数值积分法展开式处的在再作Tailorxxyyn1)()()(

!3)(2)()()(41312111hOxyhxyhxyhxyhxynnnnn)4(63211KKKhyynn因此)()(!3)(2)(4131211hOxyhxyhxyhynnnn----(8)----(9)比较

(7)、(8)两式,可知)()()(4hOheyxynnn因而三阶R-K方法（7）具有3阶精度仿真技术基础·数值积分法的系统仿真2.2数值积分法)22(643211KKKKhyynn),(111nnyxfK)2,2(1112KhyhxfKnn

),(3114hKyhxfKnn)(00xyy类似于(7)式,还可构造四阶(经典)Runge=Kutta方法)2,2(2113KhyhxfKnn)()(5hOhen-----------(10)因而方法

(10)有4阶精度仿真技术基础·数值积分法的系统仿真2.2数值积分法例1.使用高阶R-K方法计算初值问题1)0(5.002yxyy.1.0h取解:(1)使用三阶R-K方法(7)时1n201yK1210

2)21.0(KyK1025.121203))2(1.0(KKyK2555.1)4(61.032101KKKyy1111.1RK.m仿真技术基础·数值积分法的系统仿真2.2数值积分法其余结果

如下:(2)如果使用四阶R-K方法(10)时1n201yK12102)21.0(KyK1025.1nxnk1k2k3yn1.00000.10001.00001.10251.25551.11112.00000.20001.23451.37551.59451.24993.00000.300

01.56241.76372.09221.42844.00000.40002.04042.34232.86581.66645.00000.50002.77683.25874.16341.9993仿真技术基础·数值积分法的系统仿真2.2数值积分法2203)21.0(KyK1133.1230

4)1.0(KyK2351.1)22(61.0432101KKKKyy1111.1其余结果如下:nxnk1k2k3k4yn1.00000.10001.00001.10251.11331.23511.11112.00000.20

001.23461.37561.39211.56331.25003.00000.30001.56251.76391.79082.04231.42864.00000.40002.04082.34282.38

922.78051.66675.00000.50002.77773.26003.34764.00572.0000仿真技术基础·数值积分法的系统仿真2.1概述0.29650.2970.29750.2980.29850.2990.29950.30.30051.4

181.421.4221.4241.4261.428精精精精精Euler精3精RK精4精RK精仿真技术基础·数值积分法的系统仿真2.3常微分方程组和高阶微分方程的数值解法简介2.3常微分方程组和高阶微分方程的数值解法简介一、常微分方程组的数值解法下列包

含多个一阶常微分方程的初值问题称为常微分方程组的初值问题00212002212210012111)(),,,,()(),,,,()(),,,,(nnnnnnnyxyyyyxfyyxyyyyxfyyxyyyyxfy----------(

1)仿真技术基础·数值积分法的系统仿真(1)式具有n个未知函数做如下假设nyyyY21),(),(),(),(21YxfYxfYxfYxFn

002010020100)()()()(nnyyyxyxyxyxYY则(1)式化为矩阵形式00)(),(YxYYxFY----------(2)2

.3常微分方程组和高阶微分方程的数值解法简介仿真技术基础·数值积分法的系统仿真只要将以前所介绍的各种求解方法中的函数转化为函数向量,即可得到相应的常微分方程组的数值解法这里只介绍求解微分方程组的计算机实现[x,Y]=ODE45('F',xspan,Y0)首先

编制微分方程组的函数文件：functionz=F(x,Y)z=F(x,Y);然后使用求解命令ode45F为微分方程组的文件名xspan为需求值的节点向量Y0为函数向量的初值x为自变量，Y为函数值矩阵2.3常微分方程组和高阶微分方程的数值解法简介仿真技术基础·数值积分法的系统仿真例1.求解微分方程组

5.1)0(5320)0(222121211yyyxyyyyxy解:首先编制微分方程组的m文件functionz=fun(x,y)z(1)=x-y(1)+2*y(2);z(2)=2*x-3*y(1)-5*y(2)

;z=z';再编写求解程序)20(x2.3常微分方程组和高阶微分方程的数值解法简介仿真技术基础·数值积分法的系统仿真2.2数值积分法xspan=0:0.1:2;y0=[0,1.5]';[x,y]=

ode45('fun',xspan,y0)plot(x,y)00.511.52-0.500.511.5xYY1Y2运行后得仿真技术基础·数值积分法的系统仿真二、高阶常微分方程的数值解法简介例2.求下列高阶微分方程的数值解03yyyy1)0(,1)

0(,0)0(yyy)20(x解:显然yyyy3假设yy1yy2yy3则21yy32yy12333yyyy1)0(,1)0(,0)0(321yyy2.3常微分方程组和高阶微分方程的数值解法简介仿真技术基础·数值积分法的系统

仿真即高阶问题化为微分方程组的初值问题21yy32yy12333yyyy1)0(,1)0(,0)0(321yyyGaojiefangcheng.mgaojie.mfuncti

onz=gaojie(x,y)z=[y(2);y(3);y(1)*y(2)+3*y(3)];2.3常微分方程组和高阶微分方程的数值解法简介仿真技术基础·数值积分法的系统仿真xy000.10000.09450.20000.17540.30000.23810.40

000.27650.50000.28220.60000.24360.70000.14510.8000-0.03420.9000-0.32301.0000-0.7586functiongaojiefangcheng()xspan=0:0.1:1;y0=[0,

1,-1]';[x,y]=ode45('gaojie',xspan,y0);[x,y(:,1)]plot(x,y(:,1))xlabel('x')ylabel('y')2.3常微分方程组和高阶微分方程的数值解法简介仿真技术基础·数值积分法的系统

仿真00.20.40.60.81-0.8-0.6-0.4-0.200.20.40.6xy2.3常微分方程组和高阶微分方程的数值解法简介仿真技术基础·数值积分法的系统仿真2.4阿达姆斯法2.4阿达姆斯法基本思

想：在计算yn+1时，充分利用前面求出的若干点值（yn,yn-1,yn-2,yn-3,…）的信息，以期提高计算精度。设一阶常微分方程y’(t)=f[t,y(t)],其初始条件为：y(t)|t=0=y(0)=y0将方程两端从tn到tn+1求积，得：仿真技术基础·数

值积分法的系统仿真2.4阿达姆斯法根据导数f[t,y(t)]在tn,tn-1,tn-2,tn-3时的值，构造一个三次插值多项式Ｆ（t），用于替代f[t,y(t)]其中：仿真技术基础·数值积分法的系统仿真2.4阿达姆斯法具有以下特性：积分式可近似写为：仿真技术基础·数值积分法的系统仿

真2.4阿达姆斯法作变换：t=tn+xh；由于：t-tn-3=t-tn+tn–tn-3=h(x+3),t–tn-2=h(x+2),t–tn-1=h(x+1)，可以得出：仿真技术基础·数值积分法的系统仿真2.4阿达姆斯法进一步可得：四

阶阿达姆斯显式公式，也称外推公式。若用f[t,y(t)]在tn+1,tn,tn-1,tn-2时的值作插值多项式，可得：四阶阿达姆斯隐式公式，也称内推公式。仿真技术基础·数值积分法的系统仿真2.4阿达姆斯法根据插值理论，可以得出它们的截断误差分别为：显式公式：隐式公式：截断误差均为步长h的5次方幂

，具有4阶精度，与4阶龙格-库塔法相同。k阶阿达姆斯公式的一般形式为：显式：隐式：仿真技术基础·数值积分法的系统仿真2.4阿达姆斯法阿达姆斯系数表隐式显式仿真技术基础·数值积分法的系统仿真2.4阿达姆斯法阿达姆斯公式的特点：对于显式公式，每次只需计算一次右端函数；对于隐式公式，每次只需

计算二次右端函数；且均与阶数无关，因而计算效率较高。不能自启动，计算过程中改变步长困难。阿达姆斯预报-校正公式：原理：显式预报，隐式校正。4阶预报-校正公式如下：仿真技术基础·数值积分法的系统仿真2.5误差分析、收敛性和稳定性2

.5误差分析、收敛性和稳定性2.5.1误差分析整体截断误差：从初始值开始,由某个近似数值计算方法经n+1步准确计算（不含舍入误差）的解,微分方程的真解与之差。舍入误差：从初始值开始,由某个近似数值计算方法经n+1步计算，得出的实际值与该计算方法的真解之差。数值积分法

求解常微分方程的初值问题时，数值解与真解的误差由整体截断误差和舍入误差两部分组成，即仿真技术基础·数值积分法的系统仿真2.5误差分析、收敛性和稳定性2.5.2收敛性定义:一种数值积分格式,若对于任意固定的tn=t0+nh,当h→0（同时n→∞）时，由它在假定不作舍入误差的

情况下求得的yn→y(tn)，则称这种积分格式是收敛的。对于一般的单步法，数值积分格式为：定理：假定单步法（如上式）具有p阶精度，且增量函数φ(t,y,h)关于y满足李普希兹条件:非负常数又设初始值y0是准确

的，则整体截断误差：y(tn)-yn=O(hp)仿真技术基础·数值积分法的系统仿真2.5误差分析、收敛性和稳定性仿真技术基础·数值积分法的系统仿真结论:(1)整体截断误差比局部截断误差低一阶；(2)单步法是否收敛，判断其增量函数关于y是否满足李普

希兹条件。例：假定微分方程的右端函数关于y满足李普希兹条件，验证四阶龙格-库塔法收敛。2.5误差分析、收敛性和稳定性仿真技术基础·数值积分法的系统仿真2.5误差分析、收敛性和稳定性仿真技术基础·数值积分法的系统仿真2.5.2稳定性2.5误差分析、收敛性和稳定性仿真技术基础·数值积分法的系统仿真对于

模型方程欧拉公式Z变换特征方程欧拉公式的稳定条件：特征方程的根处于单位圆内，即设原微分方程的特征根=+j11(,)nnnnnnnyyhftyyyhy2.5误差分析、收敛性和稳定性仿真技术基础

·数值积分法的系统仿真稳定区域若原微分方程是稳定的，计算公式的稳定还要依赖于原微分方程的特征根和步长h。设原微分方程的时间常数为T，即步长h<2T。原微分方程是稳定的，采用欧拉公式计算是条件稳定。采用梯形公式是恒稳定的，采用RK-4是条件稳定的。

1,(0,0)T2.5误差分析、收敛性和稳定性仿真技术基础·数值积分法的系统仿真2.5误差分析、收敛性和稳定性•应用技巧：用两个显著不同的步长去计算，若所得数值结果基本相同，则该数值积分方法一般是稳定的；反之，很可能该数值积分方法不稳定。