建立系统微分方程一般步骤

1、建立系统微分方程一般步骤： （1 ）将系统划分为多个环节，确定各环节的输入及输出信号，每个环节都可考虑写一个方 程； （2）根据物理定律或通过实验等方法得出物理规律，列出各环节的原始方程式，并考虑适 当简化、线性化； （3）将各环节方程式联立，消去中间变量，最后得出只含有输入变量、输出变量以及参量 的系统方程式。 建立LRC电路的微分方程式 R i(t)匚 Ui(t) Uo(t) 用MATLAB言编程实现仿真的主要步骤是调用 MATLAB的 ODE（ Ordinary DifferentialEquation,常微分方程）解函数。MATLABt供的 常用ODE解函数如下： ? ode45 此算

2、法被推荐为首选算法； ? ode23 这是一个比ode45低阶的算法； ? ode113 用于更咼阶或大的标量计算； ? ode23t 用于解决难度适中的问题； ? ode23s 用于解决难度较大的冋题； ? ode15s 与ode23相冋，但要求的精度更咼; ? ode23tb 用于解决难度较大的冋题。 这些ODE解函数的调用格式基本相同。例如，ode45的基本调用格式为 t , x=ode45（方程函数名,tspan , x0, tol） 其中，方程函数名为描述系统状态方程的M函数名称，tspa n 般为仿真时间范围（例如， 取tspan=t0,tf, t0为起始计算时间，tf为终止计算时

3、间）;x0为系统状态变量初值， 应使该向量元素个数等于系统状态变量的个数；tol用来指定的精度，其默认值为10-3 （即 0.1%的相对误差），一般应用中可以直接采用默认值。函数返回两个结果t向量和x阵。由 于计算中采用了步长自动控制策略，因而 t向量不一定是等间隔。但是，仿真结果可以用 plot （t,x ）指令绘制出来。 例：电路如下图所示, R 1.6 ,L=2.1H , C=0.3F,初始状态是电感电流为零，电容电压 为0.2V , t=0时接入1.5V的电压，求0 t 10s时i（t） , u0（t）的值，并画出电流和电容电 压的关系曲线。 Example : Circuit is

4、shown below,R 1.6, L=2.1H, C=0.3F, Initial condition: inductanee current is 0, capacitanee voltage is 0.2V. Access 1.5V voltage when t=0, and Evaluate i(t) , u0(t) when 0 t 10s , and plot the curve between current capacita nee voltage. Ui(t) u(t) 解：(1 )根据基尔霍夫电压定律、电流定律得到系统的原始微分方程为 L 警 Ri(t) u(t) ui(t

5、) dt i(t) C du 0(t) dt (2) 消去中间变量i(t)，得电路微分方程： 2 LC d 比rcHoIH u0(t) ui(t) dtdt (3) 划高阶微分方程为一阶微分方程： 设 X1 u(t)，x2 u(t) 则有 x1 x2 1 x2(ui(t) x1 RCx2) LC (4) 编程 系统微分方程描述函数的程序：(rlcsys .m) fun cti on dx=rlcsys(t,x) dx=zeros(2,1); ui=1.5; R=1.6; L=2.1; C=0.3; dx(1)=x(2); dx(2)=1/(L*C)*(ui-x(1)-R*C*x(2); 主程序

6、( jierlc.m ) clear clc close t0=0; tfinal=10; x0=0.2 0; t,x=ode45(rlcsys,t0 tfinal,x0); C=0.3; cfu=x(:,1);%C function ifu=C*x(:,2);%i function figure(1) subplot(211) plot(t,cfu) title(C/V) xlabel(time/s) grid subplot(212) plot(t,ifu) title(L/A) xlabel(time/s) grid figure(2) plot(cfu,ifu) grid title(

7、C 、 L relation) xlabel(C) ylabel(L) 1、假设著名的 Lorenz 模型的状态方程表示为 ? x1 (t) 8x1(t ) / 3 x2(t)x3(t) ? x2(t) 10 x2 (t ) 10 x3 (t ) ? x3(t) x1(t)x2(t) 28x2(t) x3(t) 若令其初值为 x1(t) x2(t) 0， x3(t) 10 10 ，解方程。(采用 ode45() 函数) ? x1 (t) ? 8x1(t ) / 3 x2(t)x3(t) x2(t) ? 10 x2(t) 10 x3(t) x3(t) x1(t)x2(t) 28x2(t) x3(t) Lorenz model state equation is shown below: Assume that initial values are equation.(Note: ode45() 程序：( lorenz.m ) function dx=lorenz(t,x) dx=zeros(3,1); x1(t) x2(t) 0 ， x3(t) 10 10 , please solve the dx(1)=-8*x(1)/3+x(2)*x(3); dx(2)=-10*x(2)+10*x(3)