中立型 DDE
この例では、ddensd
を使用して中立型 DDE (遅延微分方程式) を解く方法を説明します。ここで、遅延は微分項に表示されます。この問題は、元々 Paul [1] により提示されました。
この方程式は次のとおりです。
.
履歴関数は、 に対し です。
この方程式は の項に時間遅延があるため、"中立型 DDE" と呼ばれます。時間遅延が の項にのみ存在する場合、方程式は時間遅延がもつ形式に応じて定数または状態依存型 DDE になります。
MATLAB® でこの方程式を解くには、遅延微分方程式ソルバー ddensd
を呼び出す前に、方程式、遅延、および履歴をコード化する必要があります。(ここで行ったように) これらをファイルの最後にローカル関数として含めることも、あるいは個別のファイルとして MATLAB パスのディレクトリに保存することもできます。
遅延のコード化
最初に、関数を記述して方程式の遅延を定義します。遅延を伴う方程式の最初の項は です。
function dy = dely(t,y) dy = t/2; end
遅延を伴う方程式の別の項は です。
function dyp = delyp(t,y) dyp = t-pi; end
この例では、 での遅延と での遅延がそれぞれ 1 つのみ存在します。より多くの遅延が存在する場合、それらの遅延を同じ関数ファイルに追加することで、関数がスカラーの代わりにベクトルを返すようにすることができます。
メモ: 関数はすべて例の最後にローカル関数として含まれます。
方程式のコード化
ここで、方程式をコード化する関数を作成します。この関数にはシグネチャ yp = ddefun(t,y,ydel,ypdel)
がなければなりません。ここでは以下のようになります。
t
は時間 (独立変数) です。y
は解 (従属変数) です。ydel
には の遅延が含まれます。ypdel
には の遅延が含まれます。
これらの入力はソルバーによって自動的に関数に渡されますが、変数名によって方程式のコード化方法が決まります。この場合、
ydel
ypdel
function yp = ddefun(t,y,ydel,ypdel) yp = 1 + y - 2*ydel^2 - ypdel; end
解の履歴のコード化
次に、解の履歴を定義する関数を作成します。解の履歴は時間 での解です。
function y = history(t) y = cos(t); end
方程式の求解
最後に、積分区間 を定義し、ddensd
ソルバーを使用して DDE を解きます。
tspan = [0 pi]; sol = ddensd(@ddefun, @dely, @delyp, @history, [0,pi]);
解のプロット
解の構造体 sol
には sol.x
フィールドと sol.y
フィールドがあり、ソルバーが受け入れる内部タイム ステップと、その時間に対応する解が含まれます ただし、deval
を使用して特定の点における解を計算することができます。
0
から pi
までにある等間隔の 20 点で解を計算します。
tn = linspace(0,pi,20); yn = deval(sol,tn);
計算された解と履歴を解析解に対してプロットします。
th = linspace(-pi,0); yh = history(th); ta = linspace(0,pi); ya = cos(ta); plot(th,yh,tn,yn,'o',ta,ya) legend('History','Numerical','Analytical','Location','NorthWest') xlabel('Time t') ylabel('Solution y') title('Example of Paul with 1 Equation and 2 Delay Functions') axis([-3.5 3.5 -1.5 1.5])
ローカル関数
ここでは、DDE ソルバー ddensd
が解を計算するために呼び出すローカル補助関数を紹介しています。あるいは、これらの関数を独自のファイルとして MATLAB パスのディレクトリに保存することもできます。
function yp = ddefun(t,y,ydel,ypdel) % equation being solved yp = 1 + y - 2*ydel^2 - ypdel; end %------------------------------------------- function dy = dely(t,y) % delay for y dy = t/2; end %------------------------------------------- function dyp = delyp(t,y) % delay for y' dyp = t-pi; end %------------------------------------------- function y = history(t) % history function for t < 0 y = cos(t); end %-------------------------------------------
参照
[1] Paul, C.A.H."A Test Set of Functional Differential Equations."Numerical Analysis Reports.No. 243.Manchester, UK:Math Department, University of Manchester, 1994.
参考
ddensd
| ddesd
| dde23
| deval