微分方程模型的定性分析与模拟2

微分方程模型的分析与模拟方法介绍二

【问题】考虑图1中的电路图。电路是由一个电容，一个电阻和一个电感器构成一个简单闭路。电路中每个元件的作用由在这个回路中的电流和电压之间的关系表示。一个理想的物理模型给出这个关系

电容：c c

dv dt

C

i =

电阻：()R R v f i = 电感：L L d i d t

L

v =

L

图1 问题2中的RLC电路图

其中 Vc 表示电容上的电压，i R 表示经过电阻的电流，L 表示电感等等。称函数f(x)为电阻的v-i 特征。在古典的RLC 电路理论中我们假设()f i i R =?，其中R 表示电阻。基尔霍夫电流律说明进入一个节点的电流值和等于流出电流之和。基尔霍夫电压率说明闭路上所有电压差之和为零。请对情形L=1，C=1和3()f i i i =-确定这个电路随时间变化的行为。

一、建立模型

这是一个典型的RLC 电路图，可以利用微分方程组来描述该系统的行为。其数学模型为：

3R R c R R c

d i d t

d v d t i i v

i ?=--??

?=?

(1) 其中c v : 电容上的电压，R i : 通过电阻的电流。

二、求解模型

1. 构造向量场（

用Mathematica 画出微分方程组的向量场）

2. 寻找平衡点

（1）求解平衡方程30

0R R c R R c

d i d t

d v d t i i v i ?=--=???==?

，可得平衡点为：（0，0）。

（2）判断该平衡点是否稳定

令3R R

c R f i i v g i ?=--??=??，可得(0,0)

(0,0)(0,0)(0,0)1110R c R

c

f f i v

g g i v A ?????????

?

-??

?== ?

???

??

?

， 方阵A

的特征值为12

λ=±，特征值的实部是正的，因此，平衡态是不稳定的。

3. 求解近似的线性系统

为了获得在平衡态附近更多的信息，我们将求解线性微分方程组

R R c c

i i A v v ??'?? ?= ? ?'????

即

R c R R c d i d t

d v d t

i v

i ?=-??

?=? (2) 利用特征值和特征向量方法求解此方程组。 对应于

12

λ=+

的特征向量为(2,1)

T

-

，因此上述方程组的一个复数解为

1()221t

e ?? ?-

??

取其实部和虚部得到两个线性无关的实值解：

22)()))t R c i u t e v t ????

?== ? ?+????

和

22)()))t R c t i v t e v ????

?== ? ?-????

。 因此，该线性系统的一般解可以写成

1212()(),,R c i c u t c v t c c R v ??

=+∈ ???

(3)

原问题的向量场

近似问题的向量场

三、数值模拟

进行数值模拟动态模型的理由就是为了获得系统行为的准确的定量信息。我们使用欧拉

方法模拟微分系统（1）。

一般地，考虑连续时间的微分方程模型

()x F x '=

其中12(,,...,)n x x x x =，12()(,,...,)n F x f f f =，带有初始条件00()x t x =。从初始条件开始，利用

()()(())x t h x t hF x t +-≈

可构造出利用欧拉方法模拟的算法如下：

算法：t(n)=n 步后的时间

1()x n =在时刻t(n)的第一状态变量 2()x n =

在时刻t(n)的第二状态变量

N=步数 T=模拟结束时间 输入：12(0),(0),(0),,t x x N T 过程：

1111222212((0))/01(1)()(()),())(1)()(()),())(1)()B egin

h T t N for n to N do B egin

x n x n hf x n x n x n x n hf x n x n t n t n h E nd E nd

←-=-+←++←++←+ 输出：

1122(1),...,()(1),...,()(1),...,()

t t N x x N x x N 欧拉方法的每一次迭代产生一个基于估计值()x t 的对于(1)x t +的估计。下面是利用更改

不同的初始条件(0)x 得到的关于电压V c 随电流变化的图形。

从上面的两个图形，我们发现，当从初始值(0)1,(0) 1.5R c i v =-=-开始，结果解曲线逆时针方向旋转地趋向原点，但是，在它接近原点时，解开始或多或少的具有周期行为，绕着原点旋转；当我们从原点附近出发，同样的行为发生了，只是现在解曲线旋转向外。这两种情形解都逼近同一条绕着原点的闭轨线。这条闭轨线称为极限环。

该系统完整的相图如下：

附录

1．模拟程序清单

" Initial Value iR=-1.1,vc=1.5,T=1000,N=10000"

t0=0;

T=1000;

Num=10000;

x0=-2.1;y0=2.50;

h=(T-t0)/Num;

n=0;

xnyn={};

For[n=0,n

Pic6=ListPlot[xnyn,PlotJoined→True,PlotStyle→{RGBColor[0,1,1]},Axes→True,AxesLabel →{"iR",""},Frame→True,AspectRatio→Automatic,PlotLabel→"iR(0)=-2.1 , vc(0)=2.5"]