实验二 MATLAB数值计算：二阶电路的时域分析
一、实验目的
在物理学和工程技术上，很多问题都可以用一个或一组常微分方程来描述，因此要解决相应的实际问题往往需要首先求解对应的微分方程(组)。在大多数情况下这些微分方程(组)通常是非线性的或者是超越方程(比如范德堡方程，波导本征值方程等)，很难解析地求解（精确解），因此往往需要使用计算机数值求解（近似解）。MATLAB作为一种强大的科学计算语言，其在数值计算和数据的可视化方面具有无以伦比的优势。在解决常微分方程(组)问题上，MATLAB就提供了多种可适用于不同场合(如刚性和非刚性问题)下的求解器(Solver)，例如ode45，ode15s，ode23，ode23s等等。本次实验将以二阶线性电路-RLC电路和二阶非线性电路-范德堡电路的时域计算为例，了解和学习使用MATLAB作为计算工具来解算复杂的微分方程，以期达到如下几个目的：
1. 熟练使用dsolve函数解析求解常微分方程；
2. 熟练运用ode45求解器数值求解常微分方程；
3. 了解状态方程的概念，能使用MATLAB对二阶电路进行计算和分析；
二、实验预备知识
1．微分方程的概念
未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程（Ordinary differential equations，简称odes）。n阶常微分方程的一般形式（隐式）为：
（1）
其中t为自变量。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，否则就是非线性微分方程，例如方程
就是非线性的。
2．常微分方程的解及MATLAB指令
一阶常微分方程与高阶微分方程可以互化，已知一个n阶常微分方程（显式）：
（2）
若令
，可将上式化为n个一阶常微分方程组：
（3）式称为状态方程，y1, y2, …,yn（即y, y¢, y¢¢, …, y(n-1) ）称为状态变量，其中y1（即y）就是常微分方程(2)式的解。(3)式中右边的函数f1、f2、…、fn代表各个状态变量的一阶导数的函数表达式，对于具体的方程它们有具体的形式，例如下列二阶非线性微分方程：
若令
，可将其改写成2个一阶微分方程组(状态方程)的形式：
因此
。
l 解析解
只有少部分的线性常微分方程可以解析地求解（即可以算出精确的解表达式），例如一阶常系数常微分方程
可以通过直接积分解出，而多数微分方程尤其是非线性方程则很难得到解析解。
有解析解的方程虽然可以手算解出，但是MATLAB也提供了dslove指令来求方程的解析解，其使用格式：
S = dsolve (‘方程1’, ‘方程2’,…,’初始条件1’，’初始条件2’…,’自变量’)
方程用字符串表示，自变量缺省值为t。1阶导数用D表示，2阶导数用D2表示，以此类推。S用于返回方程解析解的表达式。如果是求解方程组，则S为一个结构体数组，它的每个域存放方程组每一个解的表达式。
例1：求下列微分方程的解析解
>> s=dsolve(’D2y=sin(2*x)-y’,’y(0)=0’,’Dy(0)=1’,’x’);
>> simplify(s) %以最简形式显示s
ans =
-1/3*sin(x)*(-5+2*cos(x)) % 方程的解（符号表达式）
l 数值解
对于没有解析解的方程主要依靠计算机进行数值求解（得到的是近似解），例如方程就须通过计算机数值求解（结果是一系列解的数值而非表达式）。考虑n阶微分方程（2）式的数值求解，它等价于一阶常微分方程组(3)式。现将(3)式写成矩阵形式：
其中
为n个分量的列向量(Column vector)，也称状态向量
n个分量的列向量，其每个元素分别为(3)式右边的函数表达式
我们知道，微分方程要有唯一的确定解，必须给定初值条件。因此方程(4)式要有确定的解必须给定初值条件(t0为初始时刻)：
（7）
Matlab提供了ode45指令（ode是常微分方程的英文缩写）来求解方程(4)的数值解（近似解）。基本使用格式：
[tout, Yout] =ode45 (odefun, tspan, Y0,options)
其参数说明如下：
① odefun 一般是用M文件编写的函数，odefun代表函数名，由用户自己定义。函数返回值为(4)式右边的F(t,Y)=（f1, f2 ,…, fn）T。故odefun函数的返回值应是列向量，其最简单的编写格式为：
function F = odefun (t, Y) 【 作用是计算并返回4式中的F(t, Y) 】
其中 t 时间变量，为标量，代表计算进程中的某时刻点
Y 代表状态变量的列向量(即5式)
F 返回值F(t, Y)，为列向量（见6式）
ode45求解指令在计算时将会不断地在各个时间点调用odefun函数，并自动给输入参数t和Y赋值。
② tspan 指定方程的求解区间[t0, tf]，t0是初始时刻。
③ Y0 用户给定的初值条件，为n个分量的列向量，见(7)式。
④ options 可选项。一般情况下可缺省即可，若用户有特殊要求则须使用odeset指令设置options选项，具体用法可使用help odeset命令查询，此处不做要求。
⑤ tout 列向量，输出求解过程中区间[t0, tf]上各个计算点的时刻，即
, [t0, tf]上计算点的数目是由Matlab自动生成。
⑥ Yout 输出矩阵，其排列格式如下：
Yout的第1列代表的是状态变量y1在
（依次为tout的每个元素）各个时刻的值，由于y1=y，所以Yout的第一列就是待求方程的数值解，它显然是一系列离散的y(t)值：
；将Yout第一列提取出来并用plot(tout,Yout(:,1))指令即可绘制解y(t)的函数图像。Yout的第2列是y的一阶导数
的数值解，第3列是y的二阶导数
的数值解…。
三、实验内容
1. 二阶线性电路—RLC回路的零输入响应
当电路中含有二个动态元件(如电感、电容)时，建立的电路方程为二阶微分方程，这样的电路称为二阶电路，如果微分方程是线性的，则为线性电路，若为非线性方程，则是非线性电路，如范德堡电路。所谓零输入响应指的是电路中无外加的激励源，仅由动态元件初始储能所产生的响应。
考虑如图所示的RLC电路，假设电容原已充电，根据电路理论，此二阶电路的零输入响应可用如下二阶线性微分方程描述：
其中uc代表电容电压。给定初值条件：
元件参数L=0.5H, R=12.5W, 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; % 定义区间[0,1]上的时间序列
y=eval(S); % eval串演算函数，计算字串S（即9式的解）在t时刻的值
plot(t,y) % 绘制电压波形，即uc(t)的零输入响应
B. 数值解
方程(9)虽然可以解析求解，但也可以使用ode45指令来近似计算。令y1= uc, y2= duc/dt,则方程9改写成如下状态方程：
上式也可以写成(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)所示，可很明显看出电阻是非线性的。
范德堡电路的特性可由一个二阶非线性微分方程描述（参见邱关源的《电路》第5版P463-P465）：
其中
。上式就是著名的范德堡方程。该方程很难解析求解，必须借助计算机数值求解。令
，将范德堡方程改写成一阶微分方程组(状态方程)的形式：
上式也可以写成(4)式的形式。因此状态向量Y和F分别为
，
（15）
假定参数e=0.1，求解区间[0, 100]，以及初值条件
请仿照实验内容1，利用ode45指令求解范德堡方程，并绘制电流iL的波形（即电流的零输入响应图）。
四、实验任务
1. 输入实验内容1中提供的程序上机练习，熟悉dslove指令以及ode45指令的使用方法，并尝试解答思考题1。
2. 在任务1的基础上，独立编制实验内容2的计算程序，并尝试解答思考题3。
五、思考题
1. 在实验内容1中，若使用ode45指令解方程9式，如何得到电流的零输入响应波形？
2. 在实验内容1中，试研究R=1W，2W，3W，…，10W时电压Uc(t)的响应波形，并把它们绘制在一张图上。
3. 从任务2的电流波形图 (如下图) 中，你能得到什么结论？如何绘制范德堡电路的uc ~ iL函数图像（称为相图）？ （提示：
）
附：实验任务2的参考图
实验二
思考题1
函数circuit_2order_odefun.m
function F=circuit_2order_odefun(t,Y)
global L C R
F=[Y(2);-1/L/C*Y(1)-R/L*Y(2)];
脚本circuit_2order.m
global L C R
L=0.5;C=0.02;R=12.5;
Y0=[1;0];
tspan=[0,1];
[tout,Yout]=ode45('circuit_2order_odefun',tspan,Y0);
plot(tout,-C*Yout(:,2));
2. 函数circuit_2order_odefun.m
function F=circuit_2order_odefun(t,Y)
global L C R
F=[Y(2);-1/L/C*Y(1)-R/L*Y(2)];
脚本circuit.m
global L C R;
L=0.5;C=0.02;
for R=1:10
Y0=[1;0];
tspan=[0,1];
[tout,Yout]=ode45('circuit_2order_odefun',tspan,Y0);
plot(tout,Yout(:,1));
hold on;
end
3.任务二
函数VandePol_3order_odefun.m
function F=VandePol_2order_odefun(t,Y)
global a
F=[Y(2);a*(1-Y(1)*Y(1))*Y(2)-Y(1)];
脚本VandePol_3orde.m
global a;
a=0.1;
Y0=[0;1];
tspan=[0,100];
[tout,Yout]=ode45('VandePol_3order_odefun',tspan,Y0);
plot(tout,Yout（：，1）)；
u-i图像
u_i.m
global a;
a=0.1;
Y0=[0;1];
tspan=[0,100];
[tout,Yout]=ode45('VandePol_3order_odefun',tspan,Y0);
u=1/a*Yout(:,2)+1/3*Yout(:,1).^3-Yout(:,1);
plot(Yout(:,1),u)
第二篇：实验二MATLAB数值计算 (2)
实验二 MATLAB数值计算
1.实验目的
（1）掌握MATLAB变量的使用
（2）掌握MATLAB数组的创建，
（3）掌握MATLAB数组和矩阵的运算。
（4）熟悉MATLAB多项式的运用
2.实验仪器
（1）Matlab6.5应用软件安装版 一套
（3）PC机 一台
3. 实验原理
矩阵运算和数组运算在MATLAB中属于两种不同类型的运算，数组的运算是从数组元素出发，针对每个元素进行运算，矩阵的运算是从矩阵的整体出发，依照线性代数的运算规则进行。
4. 实验步骤
（1）使用冒号生成法和定数线性采样法生成一维数组。
（2）使用MATLAB提供的库函数reshape，将一维数组转换为二维和三维数组。
（3）使用逐个元素输入法生成给定变量，并对变量进行指定的算术运算、关系运算、逻辑运算。
（4）使用MATLAB绘制指定函数的曲线图，将所有输入的指令保存为M文件。
5. 实验报告内容
（1）在[0,2*pi]上产生50个等距采样数据的一维数组，用两种不同的指令实现。
a=0:2*pi/49:2*pi %b=linspace(0,2*pi,50)
a =
Columns 1 through 4
0 0.1282 0.2565 0.3847
Columns 5 through 8
0.5129 0.6411 0.7694 0.8976
Columns 9 through 12
1.0258 1.1541 1.2823 1.4105
Columns 13 through 16
1.5387 1.6670 1.7952 1.9234
Columns 17 through 20
2.0517 2.1799 2.3081 2.4363
Columns 21 through 24
2.5646 2.6928 2.8210 2.9493
Columns 25 through 28
3.0775 3.2057 3.3339 3.4622
Columns 29 through 32
3.5904 3.7186 3.8468 3.9751
Columns 33 through 36
4.1033 4.2315 4.3598 4.4880
Columns 37 through 40
4.6162 4.7444 4.8727 5.0009
Columns 41 through 44
5.1291 5.2574 5.3856 5.5138
Columns 45 through 48
5.6420 5.7703 5.8985 6.0267
Columns 49 through 50
6.1550 6.2832
>>
（2）将一维数组A=1:18，转换为2×9数组和2×3×3数组。
>> a=1:18
a =
Columns 1 through 7
1 2 3 4 5 6 7
Columns 8 through 14
8 9 10 11 12 13 14
Columns 15 through 18
15 16 17 18
>> c=reshape(1:18,2,9)
c =
Columns 1 through 7
1 3 5 7 9 11 13
2 4 6 8 10 12 14
Columns 8 through 9
15 17
16 18
>>d=reshape(1:18,2,3,3)
d(:,:,1) =
1 3 5
2 4 6
d(:,:,2) =
7 9 11
8 10 12
d(:,:,3) =
13 15 17
14 16 18
（3）A=[0 2 3 4 ;1 3 5 0],B=[1 0 5 3;1 5 0 5]，计算矩阵A转置与B乘积,数组A、B乘积，计算A&B,A|B,~A,A= =B,A>B。
a=[0 2 3 4;1 3 5 0]% M文件
b=[1 0 5 3;1 5 0 5]
d=a'*b
c=a.*b
m=(a&b)
n=(a|b)
p=(a==b)
k=(a>b)
a =
0 2 3 4
1 3 5 0
b =
1 0 5 3
1 5 0 5
d =
1 5 0 5
5 15 10 21
8 25 15 34
4 0 20 12
c =
0 0 15 12
1 15 0 0
m =
0 0 1 1
1 1 0 0
n =
1 1 1 1
1 1 1 1
p =
0 0 0 0
1 0 0 0
k =
0 1 0 1
0 0 1 0
（4）绘制y= 0.5
-t*t*sin(t),t=[0,pi]并标注峰值和峰值时间,添加标题y= 0.5-t*t*sint，将所有输入的指令保存为M文件\
a=0.5
b=1/3
t=0:0.01:pi
y=a*exp(b*t)-t.*t.*sin(t)
[y_max,t_max]=max(y)
t_text=['t=',num2str(t(t_max))]
y_text=['y=',num2str(y_max)]
max_text=char('maximum',t_text,y_text)
tit=['y=0.5*exp(',num2str(b),'t)-t*t*sin(t)']
hold on
plot(t,y,'b')
plot(t(t_max),y_max,'r.')
text(t(t_max)+0.3,y_max+0.05,max_text)
title(tit),xlabel('t'),ylabel('y'),hold off