求解刚性晶体管微分代数方程
此示例说明如何使用 ode23t
求解描述电路的刚性微分代数方程 (DAE) [1]。在示例文件 amp1dae.m 中编码的单晶体管放大器问题可以重写为半显式形式,但本示例采用它的原始形式 进行求解。该问题包括一个常奇异质量矩阵 。
晶体管放大器电路包含六个电阻器、三个电容器和一个晶体管。
初始电压信号为 。
工作电压为 。
节点处的电压由 给出。
电阻 的值是常量,通过每个电阻的电流满足 。
电容 的值是常量,通过每个电容的电流满足 。
目标是求解通过节点 5 的输出电压,。
要在 MATLAB® 中求解此方程,您需要编写方程代码,编写质量矩阵代码,并在调用求解器 ode23t
之前设置初始条件和积分区间。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
编写质量矩阵代码
根据基尔霍夫定律均衡通过每个节点(1 到 5)的电流,可以得到包含五个描述电路的方程的方程组:
此方程组的质量矩阵通过聚合方程左侧的导数项得到,其形式如下
其中,当 时,。
用适当的常量 创建一个质量矩阵,然后使用 odeset
函数指定该质量矩阵。尽管质量矩阵显然是奇异矩阵,请仍保留 'MassSingular'
选项的默认值 'maybe'
,以测试求解器对 DAE 问题的自动检测。
c = 1e-6 * (1:3);
M = zeros(5,5);
M(1,1) = -c(1);
M(1,2) = c(1);
M(2,1) = c(1);
M(2,2) = -c(1);
M(3,3) = -c(2);
M(4,4) = -c(3);
M(4,5) = c(3);
M(5,4) = c(3);
M(5,5) = -c(3);
opts = odeset('Mass',M);
编写方程代码
函数 transampdae
包含此示例的方程组。该函数定义所有电压和常量参数的值。在质量矩阵中对在方程左侧聚合的导数编码,而 transampdae
对方程的右侧编码。
function dudt = transampdae(t,u) % Define voltages and parameters Ue = @(t) 0.4*sin(200*pi*t); Ub = 6; R0 = 1000; R15 = 9000; alpha = 0.99; beta = 1e-6; Uf = 0.026; % Define system of equations f23 = beta*(exp((u(2) - u(3))/Uf) - 1); dudt = [ -(Ue(t) - u(1))/R0 -(Ub/R15 - u(2)*2/R15 - (1-alpha)*f23) -(f23 - u(3)/R15) -((Ub - u(4))/R15 - alpha*f23) (u(5)/R15) ]; end
注意:该函数作为局部函数包含在示例的末尾。
编写初始条件代码
设置初始条件。此示例对通过 [1] 中计算的每个节点的电流使用一致初始条件。
Ub = 6; u0(1) = 0; u0(2) = Ub/2; u0(3) = Ub/2; u0(4) = Ub; u0(5) = 0;
求解方程组
使用 ode23t
求出 DAE 方程组在时间区间 [0 0.05]
内的解。
tspan = [0 0.05]; [t,u] = ode23t(@transampdae,tspan,u0,opts);
绘制结果
绘制初始电压 和输出电压 。
Ue = @(t) 0.4*sin(200*pi*t); plot(t,Ue(t),'o',t,u(:,5),'.') axis([0 0.05 -3 2]); legend('Input Voltage U_e(t)','Output Voltage U_5(t)','Location','NorthWest'); title('One Transistor Amplifier DAE Problem Solved by ODE23T'); xlabel('t');
参考
[1] Hairer, E., and Gerhard Wanner.Solving Ordinary Differential Equations II:Stiff and Differential-Algebraic Problems.Springer Berlin Heidelberg, 1991, p. 377.
局部函数
此处列出 ODE 求解器 ode23t
为计算解而调用的局部辅助函数。您也可以将此函数作为单独的文件保存在 MATLAB 路径上的目录中。
function dudt = transampdae(t,u) % Define voltages and parameters Ue = @(t) 0.4*sin(200*pi*t); Ub = 6; R0 = 1000; R15 = 9000; alpha = 0.99; beta = 1e-6; Uf = 0.026; % Define system of equations f23 = beta*(exp((u(2) - u(3))/Uf) - 1); dudt = [ -(Ue(t) - u(1))/R0 -(Ub/R15 - u(2)*2/R15 - (1-alpha)*f23) -(f23 - u(3)/R15) -((Ub - u(4))/R15 - alpha*f23) (u(5)/R15) ]; end