Solve Stiff Transistor Differential Algebraic Equation

Open Live Script

This example shows how to useode23tto solve a stiff differential algebraic equation (DAE) that describes an electrical circuit [1]. The one-transistor amplifier problem coded in the example fileamp1dae.mcan be rewritten in semi-explicit form, but this example solves it in its original form $M u^{'} = ϕ (u)$ ．这一问题包括一个常数,单数马斯s matrix $M$ ．

The transistor amplifier circuit contains six resistors, three capacitors, and a transistor.

The initial voltage signal is $U_{e} (t) = 0 ． 4 \sin (200 π t)$ ．
The operating voltage is $U_{b} = 6$ ．
The voltages at the nodes are given by $U_{i} (t) (i = 1, 2, 3, 4, 5)$ ．
The values of the resistors $R_{i} (i = 1, 2, 3, 4, 5, 6)$ are constant, and the current through each resistor satisfies $I = U / R$ ．
The values of the capacitors $C_{i} (i = 1, 2, 3)$ are constant, and the current through each capacitor satisfies $I = C \cdot dU / dt$ ．

The goal is to solve for the output voltage through node 5, $U_{5} (t)$ ．

解这个方程在MATLAB®,哟u need to code the equations, code a mass matrix, and set the initial conditions and interval of integration before calling the solverode23t．You can either include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

Code Mass Matrix

Using Kirchoff's law to equalize the current through each node (1 through 5), you can obtain a system of five equations describing the circuit:

$\begin{array}{l} node 1 : \frac{U_{e} (t)}{R_{0}} - \frac{U_{1}}{R_{0}} + C_{1} ({U_{2}}^{'} - {U_{1}}^{'}) = 0, \\ node 2 : \frac{U_{b}}{R_{2}} - U_{2} (\frac{1}{R_{1}} + \frac{1}{R_{2}}) + C_{1} ({U_{1}}^{'} - {U_{2}}^{'}) - 0 ． 01 f (U_{2} - U_{3}) = 0, \\ node 3 : f (U_{2} - U_{3}) - \frac{U_{3}}{R_{3}} - C_{2} {U_{3}}^{'} = 0, \\ node 4 : \frac{U_{b}}{R_{4}} - \frac{U_{4}}{R_{4}} + C_{3} ({U_{5}}^{'} - {U_{4}}^{'}) - 0 ． 99 f (U_{2} - U_{3}) = 0, \\ node 5 : - \frac{U_{5}}{R_{5}} + C_{3} ({U_{4}}^{'} - {U_{5}}^{'}) = 0 ． \end{array}$

The mass matrix of this system, found by collecting the derivative terms on the left side of the equations, has the form

$M = (\begin{array}{ccccc} - c_{1} & c_{1} & 0 & 0 & 0 \\ c_{1} & - c_{1} & 0 & 0 & 0 \\ 0 & 0 & - c_{2} & 0 & 0 \\ 0 & 0 & 0 & - c_{3} & c_{3} \\ 0 & 0 & 0 & c_{3} & - c_{3} \end{array}),$

where $c_{k} = k \times 1 0^{- 6}$ for $k = 1, 2, 3$ ．

Create a mass matrix with the appropriate constants $c_{k}$ , and then use theodesetfunction to specify the mass matrix. Even though it is apparent that the mass matrix is singular, leave the'MassSingular'option at its default value of'maybe'to test the automatic detection of a DAE problem by the solver.

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);

Code Equations

The functiontransampdaecontains the system of equations for this example. The function defines values for all of the voltages and constant parameters. The derivatives gathered on the left side of the equations are coded in the mass matrix, andtransampdaecodes the right side of the equations.

functiondudt = transampdae(t,u)% Define voltages and parametersUe = @(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 equationsf23 = 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

Note: This function is included as a local function at the end of the example.

Code Initial Conditions

Set the initial conditions. This example uses the consistent initial conditions for the current through each node computed in [1].

Ub = 6; u0(1) = 0; u0(2) = Ub/2; u0(3) = Ub/2; u0(4) = Ub; u0(5) = 0;

Solve System of Equations

Solve the DAE system over the time interval[0 0.05]usingode23t．

tspan = [0 0.05]; [t,u] = ode23t(@transampdae,tspan,u0,opts);

Plot Results

Plot the initial voltage $U_{e} (t)$ and output voltage $U_{5} (t)$ ．

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');

Figure contains an axes object. The axes object with title One Transistor Amplifier DAE Problem Solved by ODE23T contains 2 objects of type line. These objects represent Input Voltage U_e(t), Output Voltage U_5(t).

References

[1] Hairer, E., and Gerhard Wanner.Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems．Springer Berlin Heidelberg, 1991, p. 377.

Local Functions

Listed here is the local helper function that the ODE solverode23tcalls to calculate the solution. Alternatively, you can save this function as its own file in a directory on the MATLAB path.

functiondudt = transampdae(t,u)% Define voltages and parametersUe = @(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 equationsf23 = 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