Main Content

境界層を使用した選点による非線形 ODE の解法

ライブ スクリプトを開く

この例では、Curve Fitting Toolbox™ のスプライン コマンド群を使用して、非線形常微分方程式 (ODE) を解く方法を説明します。

問題

次の非線形特異摂動の問題について考えます。

epsilon D^2g(x) + (g(x))^2 = 1 on [0..1]

Dg(0) = g(1) = 0.

この問題は epsilon = .001 でも既にかなり難しいため、控えめに次を選択します。

epsilon = .1;

近似領域

指定されたブレーク シーケンス breaks を使用する C^1 区分的 3 次式から、選点により近似解を求めます。したがって、次数 k が 4 になるようにします。

breaks = (0:4)/4;
k = 4;

対応する節点シーケンスを次のように取得します。

knots = augknt(breaks,k,2)
knots = 1×14

         0         0         0         0    0.2500    0.2500    0.5000    0.5000    0.7500    0.7500    1.0000    1.0000    1.0000    1.0000

どの次数と節点を選択しても、対応するスプライン空間の次元は次のようになります。

n = length(knots) - k
n = 10

離散化

自由度の数が 10 であれば、多項式区分ごとに 2 つのサイトで点を選ぶと予想するデータとうまく近似され、合計で 8 条件の場合は、2 つの周辺条件を追加すると条件の数が全部で 10 になります。

区間ごとに 2 つのガウス サイトを選択します。単位長さの '標準' 区間 [-1/2 .. 1/2] の場合、次の 2 つのサイトになります。

gauss = .5773502692*[-1/2; 1/2];

ここから、次を使用して選点サイトのコレクション全体を取得します。

ninterv = length(breaks)-1;
temp = (breaks(2:ninterv+1)+breaks(1:ninterv))/2;
temp = temp([1 1],:) + gauss*diff(breaks);
colsites = temp(:).';

数値的な問題

解決する数値的な問題は、任意の次数および任意の節点をもつ、非線形システムを満たす区分的多項式 (pp) y を見つけることです。

Dy(0) = 0

(y(x))^2 + epsilon D^2y(x) = 1 for x in colsites

y(1) = 0

線形化

y が解の現在の近似である場合、ニュートン法による、さらに適切と推定される (?) 解 z に対する線形問題は次のとおりです。

Dz(0) = 0

w_0(x)z(x) + epsilon D^2z(x) = b(x) for x in colsites

z(1) = 0

ここで、w_0(x) := 2y(x) および b(x) := (y(x))^2 + 1 です。

実際、w_0(1) := 1、w_1(0) := 1 および

w_2(x) := epsilon, w_1(x) := 0 for x in colsites

を選択して、他のすべての値 w_0、w_1、w_2 を選択し、b をまだゼロとして指定していない場合、システムに均一な形状を与えることができます。

w_0(x)z(x) + w_1(x)Dz(x) + w_2(x)D^2z(x) = b(x) for x in sites

ここで、

sites = [0,colsites,1];

解決する線形システム

このシステムは、その解 z の B スプライン係数のシステムに変換されます。これには、すべての sites の各 x における、すべての関連 B スプラインについての 0 階微分、1 階微分および 2 階微分が必要です。これらの値は、spcol コマンドによって提供されます。

次に、spcol のドキュメンテーションの重要部分を示します。

SPCOL B-spline collocation matrix.

COLLOC = SPCOL(KNOTS,K,TAU) is the matrix

[ D^m(i)B_j(TAU(i)) : i=1:length(TAU), j=1:length(KNOTS)-K ],

with D^m(i)B_j the m(i)-fold derivative of B_j,

B_j the j-th B-spline of order K for the knot sequence KNOTS,

TAU a sequence of sites,

both KNOTS and TAU are assumed to be nondecreasing, and

m(i) is the integer #{ j<i : TAU(j) = TAU(i) }, i.e., the

'cumulative' multiplicity of TAU(i) in TAU.

spcol を使用して次の行列を指定します。

colmat = spcol(knots,k, brk2knt(sites,3));

ここで使用する brk2knt は sites の各エントリを 3 倍にするため、sites の各 x について、colmat で、すべての関連する B スプラインの x における値、1 階微分、および 2 階微分が得られます。

ここから、重み w_0(x), w_1(x), w_2(x) を使用して x に対する 3 倍の行を組み合わせることで、線形システムの行列の x に対応する行を取得し、選点行列を取得します。

Y の初期推定の必要性

スプライン空間からの現在の近似 y も必要です。最初に、sites における pp 空間から妥当な初期推定を内挿して現在の近似を取得します。この推定には、次の放物線を使用します。

x^2 - 1

これは端点条件を満たし、またスプライン空間内に存在します。sites における内挿により、その B 型を求めます。非スパース行列 colmat から関連する内挿行列を選択します。これはいくつかの注意が必要なステップで実行されます。

intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:);
coefs = intmat\[0 colsites.*colsites-1 0].';
y = spmak(knots,coefs.');

結果をプロットし、x^2-1 そのものとなっていることを確認します。

fnplt(y,'g');
legend('Initial Guess (x^2-1)','location','NW');
axis([-0.01 1.01 -1.01 0.01]);
hold on

反復

現在の推定 y から、改善された近似解 z を得るために線形システムを作成し、解くことができました。実際、使用可能な初期推定 y を使用して、変更 z-y が指定された許容誤差より小さい場合に終了するよう反復を設定します。

tolerance = 6.e-9;

各反復における変更 z-y の最大ノルムを次の出力に示し、この Figure に反復のそれぞれも示します。

while 1
   vtau = fnval(y,colsites);
   weights = [0 1 0;
            [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)];
            1 0 0];
   colloc = zeros(n,n);
   for j = 1:n
      colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:);
   end
   coefs = colloc\[0 vtau.*vtau+1 0].';
   z = spmak(knots,coefs.');
   fnplt(z,'k');
   maxdif = max(max(abs(z.coefs-y.coefs))); 
   fprintf('maxdif = %g\n',maxdif)
   if (maxdif<tolerance), break, end
   
   % now reiterate
   y = z;
end
maxdif = 0.206695
maxdif = 0.01207
maxdif = 3.95151e-05
maxdif = 4.43216e-10

legend({'Initial Guess (x^2-1)' 'Iterates'},'location','NW');

ニュートン反復から予想されるとおり、これは 2 次収束のようになります。

イプシロンを小さくするための準備

epsilon を減らす場合は、右端点に近い境界層を作成します。これには、等間隔でないメッシュが必要です。newknt を使用して、現在の近似から適切な (細かい) メッシュを作成します。

knots = newknt(z, ninterv+1); 
breaks = knt2brk(knots);
knots = augknt(breaks,4,2);
n = length(knots)-k;

新規ブレークの選点サイト

次に、新しい breaks に対応する選点サイトを取得します。

ninterv = length(breaks)-1;
temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2);
temp = temp([1 1], :) + gauss*diff(breaks);
colsites = temp(:).';
sites = [0,colsites,1];

さらに、新しい選点行列を求めます。

colmat = spcol(knots,k, brk2knt(sites,3));

初期推定

計算した解 z への現在のスプライン空間からの内挿として、初期推定 y を求めます。得られた内挿をプロットし、現在の解に近いことを確認します。

intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:);
y = spmak(knots,[0 fnval(z,colsites) 0]/intmat.');
fnplt(y,'c');
ax = gca;
h = ax.Children;
legend(h([6 5 1]),{'Initial Guess (x^2-1)' 'Iterates' ...
                   'New Initial Guess for New Value of epsilon'}, ...
                   'location','NW');

小さいイプシロンを使用した反復

ここでは、epsilon を 3 で除算して上記の反復を繰り返します。収束は再び 2 次となります。

epsilon = epsilon/3;
while 1
   vtau = fnval(y,colsites);
   weights = [0 1 0;
            [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)];
            1 0 0];
   colloc = zeros(n,n);
   for j = 1:n
      colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:);
   end
   coefs = colloc\[0 vtau.*vtau+1 0].';
   z = spmak(knots,coefs.');
   fnplt(z,'b');
   maxdif = max(max(abs(z.coefs-y.coefs)));
   fprintf('maxdif = %g\n',maxdif)
   if (maxdif<tolerance), break, end
   
   % now reiterate
   y = z;
end
maxdif = 0.237937
maxdif = 0.0184488
maxdif = 0.000120467
maxdif = 4.78116e-09

ax = gca;
h = ax.Children;
legend(h([10 9 5 4]), ...
       {'Initial Guess (x^2-1) for epsilon = .1' 'Iterates' ...
        sprintf('Initial Guess for epsilon = %.3f',epsilon) ...
        'Iterates'}, 'location','NW');

非常に小さいイプシロン

epsilon が非常に小さい場合は、毎回 epsilon を 3 で除算してこれらの計算を繰り返すだけです。

for ee = 1:4
   knots = newknt(z, ninterv+1); 
   breaks = knt2brk(knots);
   knots = augknt(breaks,4,2); 
   n = length(knots)-k;

   ninterv = length(breaks)-1;
   temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2);
   temp = temp([1 1], :) + gauss*diff(breaks);
   colsites = temp(:).';
   sites = [0,colsites,1];

   colmat = spcol(knots,k, brk2knt(sites,3));

   intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:);
   y = spmak(knots,[0 fnval(z,colsites) 0]/intmat.');
   fnplt(y,'c')

   epsilon = epsilon/3;
   while 1
      vtau = fnval(y,colsites);
      weights = [0 1 0;
               [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)];
               1 0 0];
      colloc = zeros(n,n);
      for j = 1:n
       colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:);
      end
      coefs = colloc\[0 vtau.*vtau+1 0].';
      z = spmak(knots,coefs.');
      fnplt(z,'b');
      maxdif = max(max(abs(z.coefs-y.coefs)));
      if (maxdif<tolerance), break, end
      
      % now reiterate
      y = z;
   end
end
ax = gca;
h = ax.Children;
legend(h([30 29 25 24]), ...
       {'Initial Guess (x^2-1) for epsilon = .1' 'Iterates' ...
        'Initial Guesses for epsilon = .1/3^j, j=1:5' ...
        'Iterates'},'location','NW');

最も小さいイプシロンに使用したブレークのプロット

次に、breaks の最終的な分布を示します。この場合は newknt がうまく機能していることがわかります。

breaks = fnbrk(fn2fm(z,'pp'),'b');
bb = repmat(breaks,3,1);
cc = repmat([0;-1;NaN],1,length(breaks));
plot(bb(:),cc(:),'r');
hold off
ax = gca;
h = ax.Children;
legend(h([31 30 26 25 1]), ...
       {'Initial Guess (x^2-1) for epsilon = .1' 'Iterates' ...
        'Initial Guesses for epsilon = .1/3^j, j=1:5' ...
        'Iterates' 'Breaks for epsilon = .1/3^5'},'location','NW');

最も小さいイプシロンの残差のプロット

次の ODE を解こうとしていることを思い出してください。

epsilon D^2g(x) + (g(x))^2 = 1 on [0..1]

チェックとして、最も小さいイプシロンについて計算した解の残差を計算およびプロットします。これも満足のいく結果に見えます。

xx = linspace(0,1,201);
plot(xx, 1 - epsilon*fnval(fnder(z,2),xx) - (fnval(z,xx)).^2)
title('Residual for the Computed Solution for Smallest epsilon');