1、若方程中未知函数及其各阶导数都是一次的,称为线性常微分方程,否则就是非线性微分方程,例如方程就是非线性的。2常微分方程的解及MATLAB指令一阶常微分方程与高阶微分方程可以互化,已知一个n阶常微分方程(显式): (2)若令,可将上式化为n个一阶常微分方程组:(3)(3)式称为状态方程,y1, y2, ,yn(即y, y, y, , y(n-1) )称为状态变量,其中y1(即y)就是常微分方程(2)式的解。(3)式中右边的函数f1、f2、fn代表各个状态变量的一阶导数的函数表达式,对于具体的方程它们有具体的形式,例如下列二阶非线性微分方程:,可将其改写成2个一阶微分方程组(状态方程)的形式:因此
2、。解析解只有少部分的线性常微分方程可以解析地求解(即可以算出精确的解表达式),例如一阶常系数常微分方程可以通过直接积分解出,而多数微分方程尤其是非线性方程则很难得到解析解。有解析解的方程虽然可以手算解出,但是MATLAB也提供了dslove指令来求方程的解析解,其使用格式:S = dsolve (方程1, 方程2,初始条件1,初始条件2 ,自变量)方程用字符串表示,自变量缺省值为t。1阶导数用D表示,2阶导数用D2表示,以此类推。S用于返回方程解析解的表达式。如果是求解方程组,则S为一个结构体数组,它的每个域存放方程组每一个解的表达式。例1:求下列微分方程的解析解 s=dsolve(D2y=s
3、in(2*x)-y,y(0)=0,Dy(0)=1,x); simplify(s) %以最简形式显示sans =-1/3*sin(x)*(-5+2*cos(x) % 方程的解(符号表达式)数值解对于没有解析解的方程主要依靠计算机进行数值求解(得到的是近似解),例如方程就须通过计算机数值求解(结果是一系列解的数值而非表达式)。考虑n阶微分方程(2)式的数值求解,它等价于一阶常微分方程组(3)式。现将(3)式写成矩阵形式:其中 为n个分量的列向量(Column vector),也称状态向量 n个分量的列向量,其每个元素分别为(3)式右边的函数表达式 我们知道,微分方程要有唯一的确定解,必须给定初值条
4、件。因此方程(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
5、, 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上各个计算点的时
6、刻,即, 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回路的零输入响应当电路中含有二个动态元件(如电感、电容)时,建立的电路方程为二阶微分方程,这样的电路称为二阶电路,如果微分方
7、程是线性的,则为线性电路,若为非线性方程,则是非线性电路,如范德堡电路。所谓零输入响应指的是电路中无外加的激励源,仅由动态元件初始储能所产生的响应。 考虑如图所示的RLC电路,假设电容原已充电,根据电路理论,此二阶电路的零输入响应可用如下二阶线性微分方程描述:其中uc代表电容电压。给定初值条件:元件参数L=0.5H, R=12.5 , C=0.02F。要得到t0时的零输入响应,就必须求解(9)式。A. 解析解方程(9)可以直接求解,因此有解析解。下面改用MATLAB中的dsolve指令来求方程(9)的解析解,程序如下:S=dsolve(D2u=-R/L*Du-1/L/C*u,u(0)=1Du(
8、0)=0t); % 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*y1R/L*y2)T
9、 (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 C=0.02; R=12.5;Y0=1;0; % 列向量,初值条件tspan=0,1; % 定义
10、求解区间0,1tout,Yout=ode45(circuit_2order_odefun,tspan,Y0);plot(tout,Yout(:,1); % 绘电压波形,即uc(t)的零输入响应2.二阶非线性电路范德堡(Van de Pol)电路范德堡电路由一个线性电感、一个线性电容和一个非线性电阻构成,如图(a)所示。非线性电阻的伏安特性如图(b)所示,可很明显看出电阻是非线性的。范德堡电路的特性可由一个二阶非线性微分方程描述(参见邱关源的电路第5版P463-P465):上式就是著名的范德堡方程。该方程很难解析求解,必须借助计算机数值求解。令,将范德堡方程改写成一阶微分方程组(状态方程)的形式
11、:上式也可以写成(4)式的形式。因此状态向量Y和F分别为, (15)假定参数 =0.1,求解区间0, 100,以及初值条件请仿照实验内容1,利用ode45指令求解范德堡方程,并绘制电流iL的波形(即电流的零输入响应图)。四、实验任务1. 输入实验内容1中提供的程序上机练习,熟悉dslove指令以及ode45指令的使用方法,并尝试解答思考题1。2. 在任务1的基础上,独立编制实验内容2的计算程序,并尝试解答思考题3。五、思考题1. 在实验内容1中,若使用ode45指令解方程9式,如何得到电流的零输入响应波形?2. 在实验内容1中,试研究R=1 ,2 ,3 ,10 时电压Uc(t)的响应波形,并把
12、它们绘制在一张图上。3. 从任务2的电流波形图 (如下图) 中,你能得到什么结论?如何绘制范德堡电路的uc iL函数图像(称为相图)? (提示:)附:实验任务2的参考图实验二思考题1函数circuit_2order_odefun.mglobal L C R脚本circuit_2order.mplot(tout,-C*Yout(:,2);2. 函数circuit_2order_odefun.m脚本circuit.mglobal L C R;for R=1:10hold on;end3.任务二函数VandePol_3order_odefun.mfunction F=VandePol_2order_odefun(t,Y)global aa*(1-Y(1)*Y(1)*Y(2)-Y(1);脚本VandePol_3orde.mglobal a;a=0.1;Y0=0;1;tspan=0,100;VandePol_3order_odefunplot(tout,Yout(:,1);u-i图像u_i.mu=1/a*Yout(:,2)+1/3*Yout(:,1).3-Yout(:,1);plot(Yout(:,1),u)
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1