www.gusucode.com > matlab 案例源码 matlab代码程序 > matlab/SolveStiffDAEExample.m

    %% Solve Stiff Differential Algebraic Equation
% This example shows how to use |ode23t| to solve a stiff differential
% algebraic equation (DAE) problem that arises from an electrical circuit.
% This problem is originally from p.377 of E. Hairer and G. Wanner,
% _Solving Ordinary Differential Equations II: Stiff and
% Differential-Algebraic Problems_, 2nd ed. (Berlin: Springer, 1996).
%
% The one-transistor amplifier problem coded by <matlab:edit('amp1dae.m')
% amp1dae.m> can be rewritten in semi-explicit form, but this example
% solves it in its original form $M u' = \phi(u)$. The problem includes a
% constant, singular mass matrix that is not diagonal. The transistor
% amplifier circuit contains six resistors, three capacitors, and a
% transistor.
%
% <<../transistor_amplifier.png>>
%
% The initial voltage signal is $U_e(t)=0.4 \sin\left( 200 \pi t \right)$,
% and the other parameters are constant. The goal is to solve for the
% output voltage through node 5.
%
% Set the values of the initial and operating voltages.
Ue = @(t) 0.4*sin(200*pi*t);
Ub = 6;

%%
% Using Kirchoff's law to equalize the current through each node (1 through
% 5), you obtain a system of five equations. The mass matrix of this system
% has the form
%
% $$ M = \left( \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} \right),$$
%
% where $c_k = k \times 10^{-6}$ for $k=1,2,3$.
%
% Use the |odeset| function to pass in the mass matrix. 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);

%%
% The function |transampdae.m| contains the system of equations for this
% example. Save this function in your current folder to run the example.
%
% <include>transampdae.m</include>
%

%%
% Set the initial conditions. This example uses the consistent initial
% conditions computed by Hairer and Wanner.
u0(1) = 0;
u0(2) = Ub/2;
u0(3) = Ub/2;
u0(4) = Ub;
u0(5) = 0;

%%
% Solve the DAE system over the time interval |[0 0.05]| using |ode23t|. 
tspan = [0 0.05];
[t,u] = ode23t(@transampdae,tspan,u0,opts);

%%
% Plot the initial voltage $U_e(t)$ and output voltage $U_5(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');