3.3电路一阶微分方程的求解

带入初值y(1)可0得 Ce
一源自文库
阶
则方程的解为： yt2(et -e)
微 分
方
程
从而有： y(tn)tn 2(etn -e)
的 求 解
y1 t12(et1 e)0.345919876 y2 t22(et2 e)0.86642536
3.3
计算结果列表（ y 为n 前向欧拉法计算近似值， y(t为n )精确值）
分
for j=1:M
方 程
k1=feval(odefun,T(j),Y(j));
的
k2=feval(odefun,T(j+1),Y(j)+h*k1);
求 解
Y(j+1)=Y(j)+(h/2)*(k1+k2);
end
Function y=euler_3_3_2(t,x)
y=2/t*x+t^2*exp(t ）
阶
微 分 方 程 的
又 y0 y(1)0
求 解
tn t0nh10.1n
有： y1y0h(t20 y0t02et0)0.271828183
y2
y1
h(2 t1
y1
t12et1
)0.6
8
4
7
5
5
5
7
8
微分方程 y' 2 y是t2e一t 阶线性微分方程，
3.3
t
一阶非齐次线性微分方程
可求出其通解： yt2(et＋C)
取代一阶导数
y(tktk11) tyk(tk)y'(tk)
一 阶 微
令步长 htk，1则tk
分 方 程
y(tk 1)y(tk)y'(tk)h
的 求
解
其近似值为： yk1yky'kh
近似解的误差首先是由差商近似代替微商引起的，这种近似代替所 产生的误差称为截断误差。还有一种误差称为舍入误差，这种误差 是由于计算时数值舍入引起的。
y(tk)yk
一
0
阶
微
0.003542087 分
方
0.008327999 程
的
0.014465436 求
解
0.022061313
0.031222181
0.042054424
0.054664263
3.3
function [T Y]=Trapezia_reckon(odefun,ab,ya,M)
% odefun: 微分方程
几何意义 Euler法折线法
改进Euler法平均斜率折线法
3.3
一 阶 微 分 方 程 的 求 解
3.3
例3. 应用梯形预估-矫正法解初值问题
y'2yt2et,1t2 ,y(1 )0 t
取步长h=0.1，并把计算结果与精确解比较
一
阶
解：据前向欧拉法
y（0） k1
yk
hyk'
微 分 方
梯形预估-矫正
一
y'(t)f(y,t) tt0 y(t0)y0
阶
微 分 方
程
的
求
基本思想：
解
在初值问题存在唯一解的时间区间内，在若干个时间离
散点上，用差分方程代替微分方程，然后逐点求解差分方
程，得到各时间离散点 、t 1 …t 2处的t n函数 近似y(值t ) 、 …
y1 y2 yn
一.前向欧拉法
3.3
当两相邻离散点之间的间隔较小时，用一阶差商
yk 1yk h(tf k 1,yk 1)
在任一步长内，用一段直线
代替函数 y(的t)曲线，此直
线段的斜率等于该函数在该 步长终点的斜率。
y(t)
y0 y(t0)
0
t0
近似值
一
y3
阶
微
y2
分
y1
y(t3)
方 程
的
y(t2)
求 精确值 解
y(t1)
h
h
h
t
t1
t2
t3
3.3
注：后向欧拉法的两种处理方式
a、b：计算区间
% ya：初值 y(a)
M：等分数目
% T: 离散的时间变量
％ Y梯形公式的预估校正法解
h=(ab(2)-ab(1))/M; ％步长 T=zeros(1,M+1);Y=zeros(1,M+1);
yj1yj h 2(y'jy'(j0)1)一阶微
T=ab(1):h:ab(2); Y(1)=ya;
适合中度刚性问题的求解 求
解
若ode45失效时，可尝试
使用;
当精度较低时，计算时间 比ode15s短
当精度较低时，计算时间 比ode15s短
3.3
在用常微分方程描述一个电路的暂态过程时，往往又包含着
多个变化速度相差十分悬殊的子过程，这样一类过程就认为
具有“刚性(stiff)”，描述这类过程的微分方程称为“刚性问
(1
0.2 tk )yk
0.1
tk2etk
程 的
求 解
yk1ykh 2[f(tk,yk)f(tk1,yk(0)1)]
yk
0.0（ [5t2k yktk2etk） （tk21
y(0) k1
tk12etk1） ]
3.3
计算结果列表（ y为k 梯形预估-矫正法计算 近似值，y(tn为) 精确值）
k tk
题”。
例如，电路某一变量以e-t缓慢衰减，而另一变量以e-1000t快速 一
衰减，两变量时间常数相差很大，建立的常微分方程就具有 阶
“刚性”。
微 分
刚性问题数值解的稳定性通常被最快的模式控制，刚性问题
方 程
解答的难度就在于其快变子过程的干扰。当我们试图在慢变 的
区间上求解刚性问题时，尽管快变分量的值已衰减到微不足
0 1.0 1 1.1 2 1.2 3 1.3 4 1.4 5 1.5 6 1.6 7 1.7
yk
0 0.342377789 0.858314537 1.592749643 2.598298239 3.936444114 5.678907103 7.909209216
y(tk )
0 0.345919876 0.866642536 1.607215079 2.620359552 3.967666295 5.720961527 7.963873479
tn t0 nh10.1n
3.3
计算结果列表（y n为后向欧拉法计算近似值，
y(为t n )精确值）
n tn
0 1.0 1 1.1 2 1.2 3 1.3 4 1.4 5 1.5 6 1.6 7 1.7
yn
0 0.444282775 1.106855535 2.040960612 3.308409773 4.980911323 7.141585856 9.886697539
然后将
y
( k
0替)1 代梯形公式等式右边出现的
y k1
阶 微
校正
yk (1 )1ykh 2[f(tk,yk)f(tk 1,yk (0 )1)]
分 方 程 的
yk (n 1 1 )ykh 2[f(tk,yk)f(tk 1,yk (n 1 ))]
求 解
n0,1,2,
迭代次数
当步长h足够小，且由前向欧拉法计算的已是较好 的近似，则迭代一、二次即可
-1.013245028
-1.420624329
-1.922824060
负
三. 梯形法及其预估-矫正法
3.3
y'(t)f(t, y) y(t0)y0
t t0
改进欧拉法
用一阶差商近似地代替函数在一个步长起点和终点的
一阶导数的平均值
一 阶
微
y(tk1)h y(tk)1 2[y'(tk)y'(tk 1)]
[T Y]=Trapezia_reckon (' euler_3_3_2',[1 2],0,10)
不同求解器的特点
3.3
求解器 ode45 ode23
ode113 ode23t ode15s ode23s
求解问题
特点
非刚性
一步算法；4，5阶 Runge-Kutta算法
非刚性
一步算法；2，3阶 Runge-Kutta算法
y(tn)yn
一
0
阶
微
0.074019693 分
方
0.181886958 程
的
0.330236735 求
0.526811864 解
0.780221173
1.100143681
1.497477101
正
3.3
分析：
当步长不是很小时，前向欧拉法的精度不
一
是很高。步长取定后，步数越多，误差越
阶
微
大。
分 方
① 前向Euler法为显式，后向Euler法
为隐式——须解出yk+1.
一 阶 微
② 可用迭代法
分 方
yk+1 (n+1) = yk + hf (tk+1，yk+1(n))
程 的 求
n = 0，1，2，… 解得yk+1 ，
解
其中yk+1(0) = yk + hf (tk，yk).（结合前
向欧拉法，预报）
3.3
例2. 应用后向欧拉法解初值问题
y'2yt2et,1t2 ,y(1 )0 t
取步长h=0.1，并把计算结果与精确解比较
一
解：据后向欧拉法 yn＋1
yn
2 h(
t n1
yn1
t e 2
tn＋1
n＋1
)
阶
微 分 方
即 : yn＋1
yn
ht e 2
tn＋1
n＋1
2
1 h
程 的
求 解
t n＋1
又
y0 y(1)0
y(tn )
0 0.345919876 0.866642536 1.607215079 2.620359552 3.967666295 5.720961527 7.963873479
y(tn)yn 一
0
阶 微
-0.098362899
分 方
-0.240212999
程 的
-0.433745533
求 解
-0.688050221
程
的
求
解
二、后向欧拉法
3.3
对于给定初始条件 y(t0)y0的微分方程
y'(t)f(y(t)t,)
用一阶差商近似代替 y(t在) 一个步长终点的一阶导数，
一 阶
则原微分方程化为：
微
分
方
y(tk1)hy(tk)y'(tk1)
程 的
求 解
其近似值：
yk1yky'k1h欧拉隐式公式
3.3 后向欧拉法的几何意义：
例1. 应用前向欧拉法解初值问题
3.3
y'2yt2et,1t2,y(1)0 t
取步长h=0.1，并把计算结果与精确解比较
【思路】 用欧拉法求解常微分方程的初值问题时，首先熟练掌握欧拉公式的 一
一般形式，根据具体题目写出找出欧拉公式的迭代式，并根据初始条件和所
给步长进行迭代求解。
解：据前向欧拉法
yn＋ 1ynh(t2n yntn 2etn)
n tn
0 1.0 1 1.1 2 1.2 3 1.3 4 1.4 5 1.5 6 1.6 7 1.7
yn
0 0.271828183 0.684755578 1.276978344 2.093547688 3.187445122 4.620817846 6.466396378
y(tn )
0 0.345919876 0.866642536 1.607215079 2.620359552 3.967666295 5.720961527 7.963873479
求 解
道，但这种快速变化的干扰仍严重影响数值解的稳定性和精
度，一般地说，隐型方法比显型方法具有更大的稳定性，因
此使用隐型方法求解刚性方程组更为合适.
在MATLAB中，ode23t、ode15s、ode23s、ode23tb适合求
解刚性问题。
非刚性
多步法；变阶次的
Adams -BashforthMoulton 算法
刚性 采用梯形算法
刚性
多步法；采用了数值差 分算法
刚性
一步法；2阶Rosebrock 算法
ode23tb 刚性 隐式Runge-Kutta算法
说明 大部分场合的首选算法
使用于精度较低的情形
一 阶
微
分
计算时间比ode45短
方
程
的
分 方 程 的
梯形公式 （欧拉中点公式）
y(tk 1)y(tk)h 2[y'(tk)y'(tk 1)]
求 解
近似值：
h yk1yk2(y'ky'k1)
3.3
显然，梯形公式是隐式法，一般求 y k 需1 要解方程，
常采用迭代法，初值由显式的欧拉公式给出：
预报
yk (0 )1ykh(tfk,yk)
一
前向欧拉法的几何意义：
y(t)
yk 1ykh(tfk,yk)
在任一步长内，用一段直 线代替函数 y(t的) 曲线，此 直线段的斜率等于该函数 在该步长起点的斜率。
y0 y(t0)
3.3
y3
y2
y1
y(t3)
一
阶
微
y(t2)
分 方
y(t1)
程
的
求
h
h
h
解
0
t0
t1
t2
t3
欧拉法的几何意义：过点A0（t0，y0），A1（t1，y1），…，A n （t n，y n ），斜率分别为f（t0，y0），f（t1，y1），…，f（tn， y n）所连接的一条折线，所以欧拉法亦称为欧拉折线法。
3.3
电路暂态分析的目的是为了得到 电路的时域响应。
建立动态电路的状态方程，得到一阶微
一 阶
分方程组（或一阶微分方程），再求该
微 分
方
方程组的解。
程 的
求
解
因此暂态分析的实质就是如何获得并
且求解电路的常微分方程。
3.3
3.3 一阶微分方程的求解
一阶微分方程的求解可归结为在给定初始条件下， 求微分方程的初值问题