MATLABによる遅延微分方程式の数値計算

高安亮紀（筑波大学システム情報系）
今日の前半部分は、
の15章Numerical Solution of Delay DEsを元に作った資料です。

1. 定数の遅延をもつ遅延微分方程式

遅延微分方程式は、特に数理生物学などの応用分野で非常に頻繁に登場する。遅延微分方程式(DDE)の数値解法は常微分方程式(ODE)の初期値問題の数値解法とよく似ているが、いくつかの重要な違いもある。簡単な例で説明する。
この式はで成り立つので、もしもを知るためにはの値を知る必要があり、でを知るためにはの値が必要になる。 すなわちでをみたす履歴関数を与える必要がある。
この履歴関数は初期条件というが、（ODEでいう）初期値以上のものである。つまり区間上で定義される一般的な関数であるため、DDEの解は無限次元であり、ODEの初期値問題よりもかなり複雑である。
MATLABを用いるとdde23やddesdという関数を呼び出すだけで、数値計算ができる。まず、 ()とし、次のように呼び出す。
dde = @(t,y,Z) -0.2*Z(:,1);
lambda = -0.259171101819074;
ddehist = @(t) exp(lambda*t);
opts = ddeset('AbsTol', 1.0e-8) ;
tspan = [0, 4];
delay = 1;
sol = dde23( dde, delay, ddehist, tspan, opts )
sol = フィールドをもつ struct:
solver: 'dde23' history: @(t)exp(lambda*t) discont: [0 1 2 3] x: [0 0.3087 0.6543 1 1.4000 1.7000 2 2.4000 2.7000 3 3.4000 3.7000 4] y: [1 0.9231 0.8440 0.7717 0.6957 0.6437 0.5955 0.5369 0.4967 0.4595 0.4143 0.3833 0.3546] stats: [1×1 struct] yp: [-0.2592 -0.2392 -0.2187 -0.2000 -0.1803 -0.1668 -0.1543 -0.1391 -0.1287 -0.1191 -0.1074 -0.0993 -0.0919]
このコードは許容誤差で問題を数値的に解く。変数 y はの近似値を指し、Z はの近似値を指す。構文Z(:,1)はZの最初の列を選んでいるが、この問題ではZは実際1列しかない。そしてdeval（IVPやBVPの解を評価する関数）によって評価点で上の関数値を得ることができる。
plot(sol.x, sol.y, LineWidth=1.6)
hold on
plot(sol.x, sol.yp, '--', LineWidth=1.4)
xlabel('$t$','interpreter','latex')
ylabel('$y$ and $y''$','interpreter', 'latex')
解の精度はどの程度なのか？ここでは残差を計算し、プロットしてみる。
clf
res = []; divN = 300;
for t = linspace(0,4,divN)
[y,dy] = deval(sol,t);
if t<1
z = ddehist(t-1);
else
z = deval(sol,t-1);
end
res = [res; dy - dde(t,y,z)];
end
 
plot(linspace(0,4,divN),res,LineWidth=1.4)
xlabel('$t$','interpreter','latex')
ylabel('residual')
上の履歴関数は、パラメータλを巧みに選び、が実際にはDDEの厳密解となるようにした。特に、初期点における履歴関数と解の間が非常に滑らかである（丸め誤差がなければ、解析的である）。これは非常に珍しい例で、通常はそうはならない。なぜなら、一般にはDDEの解の右微分と、履歴関数の左微分が一致しない。実際、は存在しないかもしれないし、 履歴関数が不連続かもしれない。この初期不連続性は、高次数値計算手法にとって問題となるため、取り扱いに注意が必要であり、ODEに対する初期値問題の数値解法とDDEの数値解法のもう一つの重要な違いである。
さらに悪いことに、この不連続性は伝播する。すなわちはで不連続を持ち、3次導関数がで不連続を持つ。この単純なDDEでは、定数遅延が1つしかないので、どこに不連続点があるかを特定するのが簡単である。さらに、tが大きくなるにつれて不連続性が弱まり、微分階数が大きくなっていくので、解はどんどん滑らかになる。これは "遅れ型 " (retarded) 遅延微分方程式と呼ばれる遅延微分方程式の特徴である。
dde23は定数遅延しか扱わないので、不連続性がどこで発生するかを前もって計算することができる。複数の遅延がある場合、これは厄介な作業となる。その詳細を説明するために定数遅延が2つある場合を考えてみる。
この場合、における不連続性はで現れ、で再び現れる。そして, , ,…, の各々で現れる。
ddex1
この例はMATLABにある例で、遅延微分方程式の連立系
をで初期履歴 ()から計算している。定数遅延はベクトル[1, 0.2]で指定され、遅延微分方程式はddex1deという関数で記述され、履歴はddex1histという関数で記述されている。以下のように実行するとにおける不連続性が見られる。
sol = dde23(@ddex1de,[1, 0.2],@ddex1hist,[0, 2]);
figure;
plot(sol.x,sol.y,LineWidth=1.4)
title('An example of Wille'' and Baker.');
xlabel('time t');
ylabel('solution y');
% function s = ddex1hist(t)
% % Constant history function for DDEX1.
% s = ones(3,1);
% end
%
% function dydt = ddex1de(t,y,Z)
% % Differential equations function for DDEX1.
% ylag1 = Z(:,1);
% ylag2 = Z(:,2);
% dydt = [ ylag1(1)
% ylag1(1) + ylag2(2)
% y(2) ];
% end
さらに初期履歴ににおける不連続性だけでなく、他のより強力な不連続性が存在する場合も計算できる。すなわち最初の例で初期履歴をのような矩形波にとるとで不連続なジャンプが起こる。この初期履歴からDDEを計算すると
dde = @(t,y,Z) -0.2*Z(:,1);
ddehist = @(t) (-1).^floor(-4*t);
opts = ddeset('AbsTol',1.0e-8,'Jumps',[-1,-0.75,-0.5,-0.25,0]);
tspan = [0, 4];
delay = 1;
sol = dde23( dde, delay, ddehist, tspan, opts )
sol = フィールドをもつ struct:
solver: 'dde23' history: @(t)(-1).^floor(-4*t) discont: [-1 -0.7500 -0.5000 -0.2500 0 0.2500 0.5000 0.7500 1 1.2500 1.5000 1.7500 2 2.2500 2.5000 2.7500 3 3.2500 3.5000 3.7500 4] x: [0 0.0312 0.0625 0.1562 0.2500 0.2812 0.3125 0.4062 0.5000 0.5312 0.5625 0.6562 0.7500 0.7812 0.8125 0.9062 1 1.2500 1.5000 1.7500 2 2.2500 2.5000 2.7500 3 3.2500 3.5000 3.7500 4] y: [1 1.0035 1.0097 1.0285 1.0472 1.0438 1.0375 1.0188 1 1.0035 1.0097 1.0285 1.0472 1.0438 1.0375 1.0188 1 0.9489 0.8976 0.8465 0.7953 0.7465 0.7004 0.6568 0.6157 0.5772 0.5410 0.5071 0.4753] stats: [1×1 struct] yp: [-0.2000 0.2000 0.2000 0.2000 0.2000 -0.2000 -0.2000 -0.2000 -0.2000 0.2000 0.2000 0.2000 0.2000 -0.2000 -0.2000 -0.2000 -0.2000 -0.2094 -0.2000 -0.2094 -0.2000 -0.1898 -0.1795 -0.1693 -0.1591 -0.1493 -0.1401 -0.1314 -0.1231]
dde23には、Jumpsオプションで初期履歴の不連続点の位置を知らせることもできる。そして予測された不連続点は格納されている。
sol.discont
ans = 1×21
-1.0000 -0.7500 -0.5000 -0.2500 0 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500 3.0000 3.2500 3.5000 3.7500 4.0000
解をプロットしてみる。
fplot(ddehist,[-1,0],LineWidth=1.4)
hold on
plot(sol.x,sol.y,LineWidth=1.4)
 
plot([-1,0],[0,0],'--')
t = linspace(0,4,divN);
[~,dy] = deval(sol,t);
plot(t,dy,'--',LineWidth=1.4)
解のグラフはのグラフよりも滑らかである。この解の残差を見てみると不連続点付近で残差は悪く、解の計算が難しいことがわかる。しかし、計算を繰り返すことで残差は少しずつ小さくなっていく様子が見られる。これは遅れ型遅延微分方程式のおかげで段々と解の滑らかさが増していくことで、数値計算も安定してくるからである。
clf
divN = 50;
j_prev = 0;
for j = 0.25:0.25:4
res = [];
for t = linspace(j_prev,j,divN)
[y,dy] = deval(sol,t);
if t<1
z = ddehist(t-1);
else
z = deval(sol,t-1);
end
res = [res; dy - dde(t,y,z)];
end
semilogy(linspace(j_prev,j,divN),abs(res),LineWidth=1.4)
hold on
j_prev = j;
end

2. 状態依存の遅延をもつ遅延微分方程式

遅延微分方程式を一般的に書くと
と表され、 () を遅延関数と呼ぶ。ODEでは様々な形で非線形性が現れますが、DDEで新しいのは、遅延関数が時間依存（定数遅延はこれ）、あるいは状態依存になる可能性があることである。つまり、ではなく、と遅延関数が未知関数自体に依存することがある。このとき遅延微分方程式は状態依存の遅延をもつという。次のような単純な問題を考えてみる。
初期条件はとする。遅延関数はででは遅延はないが、で遅延が発生する。dde23は定数遅延しか扱えなかったので、この場合はddesdという関数で数値計算する。時間までの解を以下のコードでプロットする。呼び出し方はdde23と似ていて、t、y(t)、y(t/2)をzと指定し、遅延関数t/2を指定する。
% Iserles' 1994 example
% z'(t) = -0.25*z(t) + z(t/2)*(1-z(t/2)),
% z(0)=0.5, on 0 <= t <= 10ˆ3
 
dde = @(t,y,z) -0.25*y+z.*(1-z);
dely = @(t,y) t/2;
hist = @(t) 0.5;
 
tend = 1e3;
sol = ddesd( dde, dely, hist, [0,tend]);
 
t = linspace(0,tend,3*tend+1);
[y,yp] = deval(sol,t);
 
figure
plot(t,y,LineWidth=1.4)
set(gca,'fontsize',16)
xlabel('t','FontSize',16)
ylabel('y','FontSize',16)
その残差はおおむね一様に程度であることがわかる。
figure
t2 = t/2;
y2 = deval(sol,t2);
r2 = yp - (-0.25*y + y2.*(1-y2));
plot(t,r2,'k')
set(gca,'fontsize',16)
xlabel('t','FontSize',16)
ylabel('residual','FontSize',16)
この方程式はパンタグラフ方程式 (pantograph equation)といい。遅延微分方程式だが、初期履歴がいらない方程式である。
もう一つMATLABのodeexamplesのddex3 を実行してみる。
ddex3
この例では次の方程式を考えている。
遅延関数はであり、2式目に現れる。初期履歴は ()で与えられる。この遅延関数が「遅延」になっているかはこのままでは分からないが、もしも遅延関数が未知関数の未来の時刻の状態を必要とするなら、この問題は解けない。実はこの問題は解がであり、でとなるが、で遅延が消失する。
tspan = [0.1 5];
 
dde = @(t,y,Z)[y(2); -Z(2,1).*y(2)^2.*exp(1 - y(2))];
dely = @(t,y) exp(1 - y(2));
hist = @(t) [log(t); 1./t];
 
sol = ddesd(dde, dely, hist, tspan);
ta = linspace(0.1,5);
ya = hist(ta);
 
plot(ta,ya,sol.x,sol.y,'o',LineWidth=1.4)
legend('y_1 exact','y_2 exact','y_1 ddesd','y_2 ddesd')
xlabel('Time t')
ylabel('Solution y')
title('D1 Problem of Enright and Hayashi')

3. 中立型遅延微分方程式

これまでは遅延関数はyのみだったが、遅延関数がyの導関数を含む場合、すなわち
という形の場合、中立型遅延微分方程式という。この場合、不連続性は伝播するにつれて弱まることはなく、解の数値計算（と解析）はより難しくなる。次の中立型遅延微分方程式を考える。
初期履歴はとする。MATLABでこの方程式を数値計算するには、遅延微分方程式ソルバー ddensd を使うが、ソルバーを呼び出す前に、方程式、遅延、および履歴をコード化する必要がある。
tspan = [0,2.5];
 
dde = @(t,y,ydel,ypdel) y + ypdel;
dely = [];
delyp = @(t,y) t-1;
% hist = @(t) 1.0;
hist = ddensd(@(t,y,ydel,ypdel) 0, [], [], 1, [-1,0]); % ddesd( dde, dely, hist, [0,tend]);
sol = ddensd( dde, dely, delyp, hist, tspan);
 
t = linspace(0.0001,2.5,2001);
[y, yp] = deval(sol, t);
 
figure
plot(t,y,LineWidth=1.4)
set(gca,'fontsize',16)
xlabel('t','FontSize',16)
ylabel('y','FontSize',16)
[~, zp] = deval(sol,t-1);
res = yp - y - zp;
 
figure
semilogy(t,abs(res./y),LineWidth=1.4)

4. Newton法

遅延微分方程式
の特性方程式の根を数値計算するのにはニュートン法が効果的。
tol = 5e-10;
A = 0.1; B = -0.03; tau=1;
F = @(x) x - A - B*exp(-tau*x);
DF = @(x) 1 + B*tau*exp(-tau*x);
x = -10.0;
Fx = F(x);
 
display(['Before step #1, ||F||_1 = ',num2str(norm(Fx,1))])
Before step #1, ||F||_1 = 650.694
itercount = 0;
 
while (itercount<=100) && (norm(Fx,1) > tol)
DFx = DF(x);
x = x - DFx\Fx;
Fx = F(x);
display(['After step # ',num2str(itercount+1),', ||F||_1 = ',num2str(norm(Fx,1))])
itercount = itercount+1;
end
After step # 1, ||F||_1 = 237.3547 After step # 2, ||F||_1 = 85.571 After step # 3, ||F||_1 = 30.0156 After step # 4, ||F||_1 = 9.8709 After step # 5, ||F||_1 = 2.772 After step # 6, ||F||_1 = 0.51374 After step # 7, ||F||_1 = 0.031365 After step # 8, ||F||_1 = 0.00014038 After step # 9, ||F||_1 = 2.8493e-09 After step # 10, ||F||_1 = 1.7764e-15
lambda = x
lambda = -5.1683

5. スペクトル選点法（時間があれば）

まずはODEの初期値問題をスペクトル選点法で解く。解は 以下のコードはChebfunというパッケージを使う。
a = 1; x0 = 1; % Problem parameters
dom = [0 1]; % Solution interval
n = 16; % Discretisation size
t = chebpts(n, dom); % Chebyshev grid
D = diffmat(n, dom); % Differentiation matrix
I = eye(n); % Identity matrix
z = zeros(1,n-1); % Zero vector
A = D - a*I; rhs = 0*t; % Discretised operator and RHS
% Enforce boundary condition via boundary bordering:
A(1,:) = [1, z]; rhs(1) = x0;
y = A\rhs; % Solve
sol = exp(a*t)*x0; % Exact solution
err = norm(y - sol) % Error
err = 6.9829e-15
MS = 'markersize';
plot(t, sol, '-', t, y, '.', MS, 15,LineWidth=1.4)
xlabel('t','FontSize',16)
ylabel('y','FontSize',16)
状態依存の遅延をもつ遅延微分方程式（パンタグラフ方程式）
解はとなる。
a = 1; b = -8; q = 0.5; x0 = 1; % Problem parameters
tau = q*t; % Pantograph-type delay
E = barymat(tau, t); % Barycentric interpolation matrix
A = D - a*I - b*E; rhs = 0*t; % Discretised operator
A(1,:) = [1, z]; rhs(1) = x0; % Boundary conditions
y = A\rhs; % Solve
sol = @(t) -7/2*t.^3 + 21/2*t.^2 - 7*t + 1;
err = norm(y - sol(t)) % Error
err = 1.8827e-15
tt = linspace(0, 1, 1001);
plot(tt, sol(tt), '-', t, y, '.', MS, 15, LineWidth=1.4)
xlabel('t','FontSize',16)
ylabel('y','FontSize',16)
定数の遅延をもつ遅延微分方程式
% Problem parameters
n = 10; a = 1; x0 = 1; p = 0.5;
% Left discretisation Right discretisation
domL = [0, p]; domR = [p, 2*p];
tL = chebpts(n, domL); tR = chebpts(n, domR); t = [tL ; tR];
DL = diffmat(n, domL); DR = diffmat(n, domR);
I = eye(n); Z = zeros(n,n);
AL = [DL - a*I, Z]; AR = [Z, DR - a*I];
rhsL = 0*tL; rhsR = 0*tR;
% Boundary and continuity conditions:
z = zeros(1,n-1);
B = [1, z, z, 0];
C = [z, 1, -1, z];
% Assemble:
A = [B ; AL(2:n,:) ; C ; AR(2:n,:)];
rhs = [x0 ; rhsL(2:n) ; 0 ; rhsR(2:n)];
% Solve and plot:
y = A\rhs;
sol = @(t) exp(a*t)*x0;
err = norm(y - sol(t))
err = 3.1171e-14
plot(tt, sol(tt), '-', t, y, '.', MS, 15, LineWidth=1.4)
xlabel('t','FontSize',16)
ylabel('y','FontSize',16)
tauL = max(tL - p,0); tauR = tR - p; % Discrete delay
EL = barymat(tauL, tL); ER = barymat(tauR, tL); % Interpolation matrices
AL = [DL - a*EL, Z]; AR = [-ER, DR]; % Discrete operators
% Assemble:
A = [B ; AL(2:n,:) ; C ; AR(2:n,:)];
rhs = [x0 ; rhsL(2:n) ; 0 ; rhsR(2:n)];
% Solve and plot:
y = A\rhs;
sol = [a*tL + 1 ; (a*p+1)+a*(tR-p).*(.5*a*(tR-p)+1)];
err = norm(y - sol)
err = 8.5566e-15
plot(t, sol, '-', t, y, '.', MS, 15, LineWidth=1.4)
分布型の遅延をもつ遅延微分方程式
解は ().
f = @(s,t) exp((t-s));
n = 16; dom = [0 1]; a = 1/2;
t = chebpts(n, dom); I = eye(n); z = zeros(1,n-1);
D = diffmat(n, dom); % Differentiation matrix
V = f(t,t').*cumsummat(n, dom); % Volterra operator
A = D - V + a*I; rhs = 0*t; % Discrete operator and RHS
A(1,:) = [1, z]; rhs(1) = x0; % Boundary conditions
y = A\rhs;
b = sqrt(a^2-2*a+5)/2;
sol = @(t) exp(-(a+1)/2*t).*(cosh(b*t) + .5*(1-a)/b*sinh(b*t));
plot(tt, sol(tt), '-', t, y, '.', MS, 15, LineWidth=1.4), shg
xlabel('t','FontSize',16)
ylabel('y','FontSize',16)
err = norm(y - sol(t))
err = 6.6632e-15