实验二 MATLAB数值计算：二阶电路的时域分析

二、实验预备知识
1．微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为 微分方程。如果未知函数是一元函数，称为常微分方程（Ordinary differential equations，简 称 odes） 。n 阶常微分方程的一般形式（隐式）为：
F (t , y, y' , y" ,, y ( n ) ) 0
tout (t0 t1 t2
t f )' , [t0, tf]上计算点的数目是由 Matlab 自动生成。
⑥ Yout 输出矩阵，其排列格式如下：
y1 (t0 ) y (t ) 1 1 Yout y1 (t f )
y2 (t0 ) y2 (t1 ) y2 (t f )
d 2 uc dt
2

R duc L dtຫໍສະໝຸດ 1 LCuc 0
(9)
其中 uc 代表电容电压。给定初值条件：
uc t 0 1 duc i ( t ) t 0 C dt
0
t 0
(10)
元件参数 L=0.5H, R=12.5, C=0.02F。要得到 t≥0 时的零输入响应，就必须求解(9)式。 A. 解析解 方程(9)可以直接求解，因此有解析解。下面改用 MATLAB 中的 dsolve 指令来求方程(9) 的解析解，程序如下： S=dsolve('D2u=-R/L*Du-1/L/C*u','u(0)=1','Du(0)=0','t'); % S 为字符型数组（字串） ，其值为方程 9 的解表达式 L=0.5;C=0.02;R=12.5; t=0:0.01:1; y=eval(S); plot(t,y) % 元件参数
(8)
Yout 的第 1 列代表的是状态变量 y1 在 t0 , t1 , t2 ,
, t f （依次为 tout 的每个元素）各个时刻
的值，由于 y1=y，所以 Yout 的第一列就是待求方程的数值解，它显然是一系列离散的 y(t) 值：
y(ti ), i 0,1,2,... f ；将 Yout 第一列提取出来并用 plot(tout,Yout(:,1))指令即可
似解） 。基本使用格式：
（7）
Matlab 提供了 ode45 指令（ode 是常微分方程的英文缩写）来求解方程(4)的数值解（近
[tout, Yout] =ode45 (odefun, tspan, Y0,options) 其参数说明如下： ① odefun 一般是用 M 文件编写的函数，odefun 代表函数名，由用户自己定义。函数返
f1 f F 2 fn
(6)
n 个分量的列向量，其每个元素分别为(3)式右边的函数表达式
我们知道，微分方程要有唯一的确定解，必须给定初值条件。因此方程(4)式要有确定的解 必须给定初值条件(t0 为初始时刻)：
y (t0 ) y1 (t0 ) y '(t ) y (t ) 0 2 0 Y0 Y (t0 ) ( n 1) y (t0 ) yn (t0 )
实验二
一、实验目的
MATLAB 数值计算：二阶电路的时域分析
在物理学和工程技术上， 很多问题都可以用一个或一组常微分方程来描述， 因此要解决 相应的实际问题往往需要首先求解对应的微分方程(组)。 在大多数情况下这些微分方程(组) 通常是非线性的或者是超越方程(比如范德堡方程， 波导本征值方程等)， 很难解析地求解 （精 确解） ，因此往往需要使用计算机数值求解（近似解） 。MATLAB 作为一种强大的科学计算语 言， 其在数值计算和数据的可视化方面具有无以伦比的优势。 在解决常微分方程(组)问题上， MATLAB 就提供了多种可适用于不同场合(如刚性和非刚性问题)下的求解器(Solver)，例如 ode45，ode15s，ode23，ode23s 等等。本次实验将以二阶线性电路-RLC 电路和二阶非线性 电路-范德堡电路的时域计算为例， 了解和学习使用 MATLAB 作为计算工具来解算复杂的微分 方程，以期达到如下几个目的： 1. 熟练使用 dsolve 函数解析求解常微分方程； 2. 熟练运用 ode45 求解器数值求解常微分方程； 3. 了解状态方程的概念，能使用 MATLAB 对二阶电路进行计算和分析；
y '' (1 y 2 ) y ' y 0
若令 y1 y, y2 y ' ，可将其改写成 2 个一阶微分方程组(状态方程)的形式：
y1 ' y2 2 y2 ' (1 y1 ) y2 y1
因此 f1 y2 ; 解析解 只有少部分的线性常微分方程可以解析地求解（即可以算出精确的解表达式） ，例如一 阶常系数常微分方程 dy dt y 1 可以通过直接积分解出，而多数微分方程尤其是非线性 方程则很难得到解析解。 有解析解的方程虽然可以手算解出， 但是 MATLAB 也提供了 dslove 指令来求方程的解析 解，其使用格式： S = dsolve (‘方程 1’, ‘方程 2’,…,’初始条件 1’，’初始条件 2’ …,’自变量’) 方程用字符串表示， 自变量缺省值为 t。 1 阶导数用 D 表示， 2 阶导数用 D2 表示， 以此类推。 S 用于返回方程解析解的表达式。如果是求解方程组，则 S 为一个结构体数组，它的每个域 存放方程组每一个解的表达式。 例 1：求下列微分方程的解析解
Y '(t ) F (t, Y (t ))
其中
(4)
y1 (t ) y y (t ) y ' Y (t ) 2 = ( n 1) yn (t ) y
(5)
为 n 个分量的列向量(Column vector)，也称状态变量
y ( n ) f (t , y ' , y" ,, y ( n 1) )
若令 y1 y, y2 y ', y3 y '',...., yn y
( n 1)
（2）
，可将上式化为 n 个一阶常微分方程组：
' y1 f1 (t , y1 , y2 ,... yn ) ' y2 f 2 (t , y1 , y2 ,... yn ) y ' f (t , y , y ,... y ) n n 1 2 n
y1 ' y2 1 R y2 ' y1 y2 LC L
(11)
上式也可以写成(4)式的形式，因此状态向量 Y=（y1，y2）T， F =（f1，f2）T =（y2，﹣1/L/C*y1﹣R/L*y2）T （12）
第一步：通过 M 文件创建 ode 函数，函数名 circuit_2order_odefun。 function F=circuit_2order_odefun(t,Y) global L C R % 定义全局变量 L、C、R，以实现参数在 MATLAB 的基本工作空% 间和函数的专用空间之间数据的传递 F=[Y(2);-1/L/C*Y(1)-R/L*Y(2)]; % 函数返回值 F, 列向量，见（12）式 注意，函数创建完后，必须存盘，存储文件名和函数名应一致！ 第二步：创建 M 脚本文件 circuit_2order.m L=0.5; C=0.02; R=12.5; % 元件参数 Y0=[1;0]; % 列向量，初值条件 tspan=[0,1]; % 定义求解区间[0,1] [tout,Yout]=ode45('circuit_2order_odefun',tspan,Y0); plot(tout,Yout(:,1)); % 绘电压波形，即 uc(t)的零输入响应 2. 二阶非线性电路—范德堡(Van de Pol)电路 范德堡电路由一个线性电感、一个线性电容和一个非线性电阻构成，如图(a)所示。非 线性电阻的伏安特性如图(b)所示，可很明显看出电阻是非线性的。
回值为(4)式右边的 F(t,Y)=（f1, f2 ,…, fn）T。故 odefun 函数的返回值应是列向量，其最简单 的编写格式为： function F = odefun (t, Y) 其中 t 【 作用是计算并返回 4 式中的 F(t, Y) 】
时间变量，为标量，代表计算进程中的某时刻点
Y 代表状态变量的列向量(即 5 式) F 返回值 F(t, Y)，为列向量（见 6 式） ode45 求解指令在计算时将会不断地在各个时间点调用 odefun 函数， 并自动给输入参数 t 和 Y 赋值。 ② tspan ③ Y0 ④ options 指定方程的求解区间[t0, tf]，t0 是初始时刻。 用户给定的初值条件，为 n 个分量的列向量，见(7)式。 可选项。一般情况下可缺省即可，若用户有特殊要求则须使用 odeset 指令设 置 options 选项，具体用法可使用 help odeset 命令查询，此处不做要求。 ⑤ tout 列向量，输出求解过程中区间[t0, tf]上各个计算点的时刻，即
（ 3）
（3）式称为状态方程，y1, y2, …,yn（即 y, y, y, …, y(n-1) ）称为状态变量，其中 y1（即 y） 就是常微分方程(2)式的解。(3)式中右边的函数 f1、f2、…、fn 代表各个状态变量的一阶导
数的函数表达式，对于具体的方程它们有具体的形式，例如下列二阶非线性微分方程：
yn (t0 ) y (t0 ) yn (t1 ) = y (t1 ) yn (t f ) y (t ) f
y '(t0 ) y '(t1 ) y '(t f )
y ( n 1) (t0 ) y ( n 1) (t1 ) ( n 1) y (t f )