Main Content

トランジスタのスティッフな微分代数方程式の求解

この例では ode23t を使用して電気回路を記述するスティッフな微分代数方程式 (DAE) を解く方法を説明します [1]。例のファイル amp1dae.m にコード記述されている 1 石トランジスタ アンプの問題は、半陽的な形式に書き換えることもできますが、この例では元の形式 Mu=ϕ(u) で解きます。この問題には、定数の特異質量行列 M が含まれます。

このトランジスタ アンプ回路には 6 つの抵抗、3 つのコンデンサ、1 つのトランジスタが含まれます。

  • 初期の電圧信号は Ue(t)=0.4sin(200πt) です。

  • 動作電圧は Ub=6 です。

  • ノードの電圧は Ui(t) (i=1,2,3,4,5) によって与えられます。

  • 抵抗器の値 Ri (i=1,2,3,4,5,6) は定数であり、各抵抗器を通る電流は I=U/R を満たします。

  • コンデンサの値 Ci (i=1,2,3) は定数であり、各コンデンサを通る電流は I=CdU/dt を満たします。

目的はノード 5 にかかる出力電圧 U5(t) を求めることです。

MATLAB® でこの方程式を解くには、方程式をコード化し、質量行列をコード化し、初期条件と積分区間を設定してからソルバー ode23t を呼び出します。(ここで行ったように) 必要な関数をファイルの最後にローカル関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パスのディレクトリに保存することもできます。

質量行列のコード化

キルヒホッフの法則を使用して各ノード (1 ~ 5) を通る電流を等しくすると、回路を記述する 5 つの方程式からなる系が得られます。

node1:  Ue(t)R0-U1R0+C1(U2-U1)=0,node2:  UbR2-U2(1R1+1R2)+C1(U1-U2)-0.01f(U2-U3)=0,node3:f(U2-U3)-U3R3-C2U3=0,node4:  UbR4-U4R4+C3(U5-U4)-0.99f(U2-U3)=0,node5:-U5R5+C3(U4-U5)=0.

この方程式系の質量行列は、方程式の左辺の微分項を集めることによって求められ、次の形式になります。

M=(-c1c1000c1-c100000-c200000-c3c3000c3-c3),

ここで、k=1,2,3 に対して ck=k×10-6 です。

適切な定数 ck を使って質量行列を作成し、関数 odeset を使用して質量行列を指定します。質量行列が特異行列であることは明らかですが、ソルバーによる DAE 問題の自動検出をテストするために、'MassSingular' オプションを既定値 'maybe' のままにします。

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 を使用して、時間区間 [0 0.05] での DAE 系の解を求めます。

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

結果のプロット

初期電圧 Ue(t) と出力電圧 U5(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, xlabel t contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Input Voltage U_e(t), Output Voltage U_5(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

参考

|

関連するトピック