计算机仿真技术-Chapter2

13
数值积分解法
数值解法基本思路


高阶微分方程，可以转换为一个一阶微分 方程组，从而归结为一阶微分方程的数值 求解 对于 dy f ( t , y ) dt y ( t 0 ) y0 求解基本思想：离散化和递推化
14
数值积分解法
微分方程的离散化
假设y在各个时间节点
a t0 t1 t2 t N 1 t N b
计算机仿真技术
Computer Simulation
1
第二章电气系统建模与常微分 方程解法
一、电气系统建模
二、数值积分算法
三、Matlab求解函数
四、Matlab绘图
2
一、电气系统建模
3
电气系统建模
电气系统的常用模型
微分方程模型
传递函数模型
状态方程模型
4
电气系统建模
电气系统的常用模型
28
三、Matlab求解函数
29
Matlab求解函数
微分方程的解析解

‘L*C*D2y+(L/R0+R*C)*Dy+(R/R0+1)y=udc‘
dsolve('eq1','eq2',...,'cond1','cond2',...,'v')
R = subs(S, new)
30
Matlab求解函数
d 2u L du R LC ( RC ) ( 1)u udc 2 dt R0 dt R0
6
电气系统建模
微分方程模型
d 2u L du R LC ( RC ) ( 1)u udc 2 dt R0 dt R0
在Matlab中建立该模型： 建立的是符号模型，可用于求解析解 以Dny…D2y,Dy表示微分 以’…’括住微分方程式 'L*C*D2y+(L/R0+R*C)*Dy+(R/R0+1)y=udc'
22
数值积分解法
龙格－库塔(Runge-Kutta)法
龙格－库塔法的思想： 用几个点上的一阶导数的线性组合来近似代替某 一点上的各阶导数 用泰勒级数展开式确定线性组合中的各加权系数。 在h内，用几个点的加权平均值来求h内的面积！ 龙格－库塔法的特点： 无需求导，递推公式
23
数值积分解法
欧拉法（矩形面积）
假定时间节点等分，即步长为h=(b-a)/N

tk 1
tk
f (t , y)dt hf (tk , yk )
算法简单，计算量小 从初值开始递推，仅一步递推就可求出后一点值， 属于单步法，可自启动。 具有一阶精度，精度较差。 是最经典的算法，是其它方法的基础。
34
Matlab求解函数
微分方程的数值解
求解算法


ode45：四/五阶龙格－库塔法，适用于大多数连续 或离散系统，但不适用于刚性系统。 ode23：二/三阶龙格－库塔法，它在误差限要求不 高和求解的问题不太难的情况下，可能会比ode45 更有效。 ode113：一种阶数可变的解法，它在误差容许要求 严格的情况下通常比ode45有效。ode113是一种多 步解法。
状态方程模型
di Ri L u udc dt du u iC dt R0
di 1 (udc Ri u ) i dt L du 1 u u (i ) dt C R0
10
电气系统建模
状态方程模型
di dt i X , X , X AX BU , y CX DU du u dt 1 R 1 L L , B L , U udc , C [0 1], D 0 A 1 1 0 C CR0 y u CX
38
四、Matlab绘图
39
Matlab绘图
用到的几个命令
clear clc hold （on/off） grid（on/off）
舍入误差：由计算机字长决定。
累积误差：以上两项误差随计算过程的进行所造
成的累积。
P阶精度
截断误差与步长h有关，若某方法的截断误差与hp+1 成正比，则称之为p阶精度。
21
数值积分解法
精度、步长
欧拉法具有一阶精度。
要想获得高精度，h需要较小
步长与精度的关系
截断误差要求h要小，舍入误差要求h要大。若h过小， 则计算步数增加，累积次数增多，累积误差就增大。 因此形成矛盾，折衷选取适当步长。
8
电气系统建模
传递函数模型
U ( s) 1 U dc ( s) LCs 2 ( L RC ) s ( R 1) R0 R0
在Matlab中建立该模型： 利用tf(num,den)函数 num=[1]，den=[L*C,L/R0+R*C,R/R0+1]
9
电气系统建模
h y (t k 1 ) y k 1 y k ( K1 2 K 2 2 K 3 K 4 ) 6 K1 f (t k , y k ) h h K 2 f (t k , y k K1 ) 2 2 h h K 3 f (t k , y k K 2 ) 2 2 K 4 f (t k h, y k hK3 )
35
Matlab求解函数
微分方程的数值解
求解算法


ode15s：一种基于数字微分公式的解法（NDFs）。 多步解法。适用于刚性系统，当用户估计要解决的 问题是比较困难的，或者不能使用ode45，或者即 使使用效果也不好，就可以用ode15s。 ode23s：半隐式龙格－库塔法，专门应用于刚性系 统，在弱误差允许下的效果好于ode15s。它能解决 某些ode15s所不能有效解决的刚性系统问题。
传递函数的响应


num=[1]，den=[L*C,L/R0+R*C,R/R0+1] sys=tf(num,den) lsim(sys,u,t)
31
Matlab求解函数
微分方程的数值解 [T,Y] = solver(‘odefun’,ts,y 0,options)
自变 量值 函数 值 ode45 ode23 ode113 ode15s ode23s ode23t ode23tb [t0,tf] 由待解 方程写 t0、 tf为自 成的m函 变量的初 值和终值 数文件 名 函数 的初 值
i R udc L C + u 用途：滤波器等 边界条件：
R0
输入：电源电压udc，负载R0 状态(输出)：电容电压u，电感电流i
5
电气系统建模
微分方程模型
di Ri L u udc dt
du u iC dt R0
整理后：
L d du u du u (C ) R(C ) u udc dt dt R0 dt R0
36
Matlab求解函数
微分方程的数值解
求解算法


ode23t：梯形积分法的一种自由插值实现。这种解 法适用于求解适度刚性的问题而用户又需要一个无 数字振荡的解法的情况。 ode23tb：梯形法与Gear法的结合
37
Matlab求解函数
微分方程的数值解

示例
di 1 (udc Ri u ) i dt L du 1 u u (i ) dt C R0
龙格－库塔(Runge-Kutta)法
二阶R－K公式
y (tk 1 ) yk 1 h yk ( K1 K 2 ) 2
K1 f (tk , yk ) K 2 f (tk h, yk K1h)
24
数值积分解法
龙格－库塔(Runge-Kutta)法
四阶R－K公式（更为常用）
上对应的数值为
y (t 0 ), y(t1 ), y (t 2 ),, y (t N 1 ), y (t N )
15
数值积分解法
微分方程的递推化
y (t1 ) y (t0 )
t1
t0 t2
f (t , y )dt f (t , y )dt
y (t2 ) y (t1 )
dy f (t , y ) dt
33
Matlab求解函数
微分方程的数值解

对于高阶微分方程，需要转换成一阶微分方程组的 形式 此时dy、y分别为n×1的向量，其中的分量分别表 示为dy(k)、y(k) function dy=odefun(t,y) dy=zeros(n,1); dy(1)=f1(y(1),··,y(n)); ·· ·· dyk ┇ f k (t , y) dt dy(n)=fn(y(1), ··,y(n)); ·· ··
26
数值积分解法
多步算法
常用线性多步法，由yn-r至yn和fn-r至fn的一次项构造 线性公式
y (tn 1 ) yn 1 i yn i h i f n i
i 0 i 0 r r
27
数值积分解法
用于病态系统的吉尔(Gear)法
前述方法都是基于固定步长，步长需限定在系统最 小时间常数的数量级 对于病态系统（如电力电子电路），系统的过渡过 程决定于最大时间常数 前述方法计算量大、误差大、稳定性差 Gear提出的Gear法对于病态系统求解比较有效
7
电气系统建模
传递函数模型
d 2u L du R LC ( RC ) ( 1)u udc 2 dt R0 dt R0
拉氏变换后：
L R LCs U (s) ( RC ) sU (s) ( 1)U (s) U dc ( s) R0 R0