Main Content

非線形 ODE の解法

このセクションでは、非線形 ODE 問題の次の点について説明します。

この例を実行できます。「Solving a Nonlinear ODE with a Boundary Layer by Collocation」。

問題

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

εD2g(x)+(g(x))2=1on[0..1]

Dg(0)=g(1)=0

近似領域

適切なブレーク シーケンスを使用する C1 区分的 3 次式から選点により近似解を求めます。たとえば、次のようにします。

breaks = (0:4)/4;

3 次式は次数 4 であるため、以下のようになります。

k = 4;

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

knots = augknt(breaks,k,2);

これによって、0 と 1 の両方で 4 倍の節点が与えられます。これは 3 次式、つまり次数 4 があることと一致しています。

これは、次があることを意味します。

n = length(knots)-k;
n = 10;

すなわち、10 の自由度です。

離散化

多項式区分ごとに 2 つのサイト (全部で 8 つのサイト) で点を選びます。これは 2 つの周辺条件と併せて、自由度 10 と一致する 10 の条件を与えます。

区間ごとに 2 つのガウス サイトを選択します。長さ 1 の "標準" 間隔 [–0.5,0.5] の場合、次の 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(:).';

数値的な問題

これを使用して、解決する数値的な問題は、非線形システムを満たす yS4,knots を見つけることです。

Dy(0)=0(y(x))2+εD2y(x)=1 for x  colsitesy(1)=0

線形化

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

Dz(0)=0w0(x)z(x)+εD2z(x)=b(x) for x  colsitesz(1)=0

ここでは、w0(x)=2y(x),b(x)=(y(x))2+1 です。実際、次を選択します。

w0(1):=1, w1(0):=1w1(x):=0, w2(x):=ε for x colsites

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

w0(x)z(x)+w1(x)Dz(x)+w2(x)D2z(x)=b(x),forx  sites

に対し

sites = [0,colsites,1]; 

解決する線形システム

z∊S4,knots であるため、この最後のシステムを z の B スプライン係数のシステムに変換します。これには、すべての x∊ "サイト" とすべての関連 B スプラインについての値、1 階微分および 2 階微分が必要です。コマンド spcol は、この目的のために明示的に記述されています。

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

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

ここから、重み w0(x),w1(x),w2(x) を使用して x に対する colmat の 3 倍の行を組み合わせることで実際の行列の x に対する行を取得し、選点行列を取得します。このためには、現在の近似 y が必要です。最初に、sites における区分的多項式空間から妥当な初期推定を内挿して現在の近似を取得します。初期推定として端点条件を満たす放物線 x2–1 を使用し、非スパース行列 colmat から行列を選択します。これはいくつかの注意が必要なステップで実行されます。

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

初期推定をプロットし、以降のプロット用にホールドをオンにします。

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 が "十分に小さい" 場合に終了するよう反復を設定します。比較的中程度の ε = .1 を選択します。

tolerance = 6.e-9;
epsilon = .1;
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
legend({'Initial Guess (x^2-1)' 'Iterates'},'location','NW');

誤差の結果の表示は次のようになります。

maxdif = 0.206695
maxdif = 0.01207
maxdif = 3.95151e-005
maxdif = 4.43216e-010

ε を減らす場合は、右端点に近い境界層を作成します。これには等間隔でないメッシュが必要です。

newknt を使用して、現在の近似から適切な細かいメッシュを作成します。

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); 
colpnts = temp(:).'; 
sites = [0,colpnts,1];

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

colmat = spcol(knots,k,sort([sites sites sites]));

現在の近似解 z を初期推定として使用します。

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

このように設定し、ε を 3 で除算して、次のステートメントから始まる以前の計算を繰り返します。

tolerance=1.e-9; 
while 1 
  vtau=fnval(y,colpnts);
  .
  .
  .

このプロセスを繰り返し実行すると、ε = 1/10、1/30、1/90、1/270、1/810 に対する解のシーケンスが生成されます。結果の解は、プロット例に示すように、0 でよりフラットになり、1 でより急になります。このプロットには、最終的なブレーク シーケンスも縦棒のシーケンスとして示されています。プロットを表示するには、「Solving a Nonlinear ODE with a Boundary Layer by Collocation」の例を実行します。

この例では、少なくとも newknt の動作は問題ありません。