このページの翻訳は最新ではありません。ここをクリックして,英語の最新版を参照してください。
境界層を使用した選点による非線形颂歌の解法
この例では,曲线拟合工具箱™のスプラインコマンド群を使用して,非線形常微分方程式(ODE)を解く方法を説明します。
問題
次の非線形特異摂動の問題について考えます。
在[0..1]上D^2g(x) + (g(x))^2 = 1
Dg(0) = g(1) = 0。
この問題はε=措施でも既にかなり難しいため,控えめに次を選択します。
ε= 1;
近似領域
指定されたブレークシーケンス休息时间を使用するC ^ 1区分的3次式から,選点により近似解を求めます。したがって,次数kが4になるようにします。
休息时间= (0:4)/ 4;k = 4;
対応する節点シーケンスを次のように取得します。
节= augknt(休息,k, 2)
节=1×14000 0.2500 0.2500 0.5000 0.5000 0.7500 0.7500 1.0000 1.0000 1.0000 1.0000
どの次数と節点を選択しても,対応するスプライン空間の次元は次のようになります。
N =长度(节)- k
n = 10
離散化
自由度の数が10であれば,多項式区分ごとに2つのサイトで点を選ぶと予想するデータとうまく近似され,合計で8条件の場合は,2つの周辺条件を追加すると条件の数が全部で10になります。
区間ごとに2つのガウスサイトを選択します。単位長さの”標準”区間[1/2 . .2 1/2)の場合,次のつのサイトになります。
高斯= .5773502692 * [1/2;1/2);
ここから,次を使用して選点サイトのコレクション全体を取得します。
ninterv =长度(减免)1;temp =(优惠(2:ninterv + 1) +优惠(1:ninterv)) / 2;Temp = Temp ([1 1]，:) + gauss*diff(breaks); / /去掉所有的空colsites =临时(:)。';
数値的な問題
解決する数値的な問題は,任意の次数および任意の節点をもつ,非線形システムを満たす区分的多項式(页）yを見つけることです。
Dy (0) = 0
(y(x))^2 + D^2y(x) = 1
y (1) = 0
線形化
yが解の現在の近似である場合,ニュートン法による,さらに適切と推定される(?)解zに対する線形問題は次のとおりです。
Dz (0) = 0
w_0(x)z(x) + D^2z(x) = b(x)对于colsites中的x
z (1) = 0
ここで,w_0 (x): y = 2 (x)およびB (x) = (y(x))^2 + 1です。
実際、w_0 (1): = 1、w_1 (0): = 1および
W_2 (x):=， w_1(x):= 0
を選択して,他のすべての値w_0、w_1、w_2を選択し,bをまだゼロとして指定していない場合,システムに均一な形状を与えることができます。
w_0 (x) z (x) + w_1 (x) Dz (x) + w_2 (x) ^ 2 z D b (x) = x (x)的网站
ここで,
colsites网站= [0,1];
解決する線形システム
このシステムは,その解zのBスプライン係数のシステムに変換されます。これには,すべての网站の各xにおける,すべての関連Bスプラインについての0階微分,1階微分および2階微分が必要です。これらの値は,spcolコマンドによって提供されます。
次に,spcolのドキュメンテーションの重要部分を示します。
SPCOL b样条配置矩阵。
COLLOC = SPCOL(KNOTS,K,TAU)为矩阵
[D ^ m (i) B_j(τ(i)): i = 1:长度(τ),j = 1:长度(节)- k),
D^m(i)B_j B_j的m(i)倍导数，
B_j结点序列的第j条b样条的K阶
TAU的一系列地点，
假设节和TAU都是非递减的，并且
m(i)为整数#{j
TAU(i)的“累积”多重性。
spcolを使用して次の行列を指定します。
Colmat = spcol(结点，k, brk2knt(sites,3));
ここで使用するbrk2kntは网站の各エントリを3倍にするため,网站の各xについて,colmat,ですべての関連するBスプラインのxにおける値,1階微分,および2階微分が得られます。
ここから,重みw_0 (x) w_1 (x) w_2 (x)を使用してxに対する3倍の行を組み合わせることで,線形システムの行列のxに対応する行を取得し,選点行列を取得します。
Yの初期推定の必要性
スプライン空間からの現在の近似yも必要です。最初に,网站における页空間から妥当な初期推定を内挿して現在の近似を取得します。この推定には,次の放物線を使用します。
x ^ 2 - 1
これは端点条件を満たし,またスプライン空間内に存在します。网站における内挿により,そのB型を求めます。非スパース行列colmatから関連する内挿行列を選択します。これはいくつかの注意が必要なステップで実行されます。
Intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3]，:);Coefs = intmat\[0 colsites。* colsites-1 0]。”;y = spmak(节,系数。');
結果をプロットし,x ^ 2 - 1そのものとなっていることを確認します。
fnplt (y,‘g’）;传奇(“初始猜测(x ^ 2 - 1)”，“位置”，“西北”）;轴([-0.01 1.01 -1.01 0.01]);持有在
反復
現在の推定yから,改善された近似解zを得るために線形システムを作成し,解くことができました。実際,使用可能な初期推定yを使用して,変更z-yが指定された許容誤差より小さい場合に終了するよう反復を設定します。
公差= 6. e-9;
各反復における変更z-yの最大ノルムを次の出力に示し,この图に反復のそれぞれも示します。
而1 vtau = fnval(y,colsites);Weights = [0 1 0;(2 * vtau。0 (n - 1) repmat(ε- 1));1 0 0];colloc = 0 (n, n);为j = 1: n colloc (j:) =重量(j:) * colmat (3 * (j - 1) + (1:3):);结束Coefs = colloc\[0 vtau.]* vtau + 1 0]。';z = spmak(节,系数。');fnplt (z,“k”）;maxdif = max (max (abs (z.coefs-y.coefs)));流(“maxdif = % g \ n”maxdif)如果(maxdif <公差),打破，结束%现在重申y = z;结束
Maxdif = 0.206695 Maxdif = 0.01207 Maxdif = 3.95151e-05 Maxdif = 4.43216e-10
传奇({“初始猜测(x ^ 2 - 1)”“迭代”}，“位置”，“西北”）;
ニュートン反復から予想されるとおり,これは2次収束のようになります。
イプシロンを小さくするための準備
εを減らす場合は,右端点に近い境界層を作成します。これには,等間隔でないメッシュが必要です。newkntを使用して,現在の近似から適切な(細かい)メッシュを作成します。
结点= newknt(z, ninterv+1);休息= knt2brk(节);节= augknt(休息,4,2);n =长度(节)- k;
新規ブレークの選点サイト
次に,新しい休息时间に対応する選点サイトを取得します。
ninterv =长度(减免)1;temp =((优惠(2:ninterv + 1) +优惠(1:ninterv)) / 2);Temp = Temp ([1 1]，:) + gauss*diff(breaks); / /去掉所有的空colsites =临时(:)。';colsites网站= [0,1];
さらに,新しい選点行列を求めます。
Colmat = spcol(结点，k, brk2knt(sites,3));
初期推定
計算した解zへの現在のスプライン空間からの内挿として,初期推定yを求めます。得られた内挿をプロットし,現在の解に近いことを確認します。
Intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3]，:);Y = spmak(结点，[0 fnval(z,colsites) 0]/intmat.');fnplt (y,“c”）;甘氨胆酸ax =;h = ax.Children;传奇(h ([6 5 1]), {“初始猜测(x ^ 2 - 1)”“迭代”．．．对新的值的新的初始猜测}，．．．“位置”，“西北”）;
小さいイプシロンを使用した反復
ここでは,εを3で除算して上記の反復を繰り返します。収束は再び2次となります。
ε=ε/ 3;而1 vtau = fnval(y,colsites);Weights = [0 1 0;(2 * vtau。0 (n - 1) repmat(ε- 1));1 0 0];colloc = 0 (n, n);为j = 1: n colloc (j:) =重量(j:) * colmat (3 * (j - 1) + (1:3):);结束Coefs = colloc\[0 vtau.]* vtau + 1 0]。';z = spmak(节,系数。');fnplt (z,“b”）;maxdif = max (max (abs (z.coefs-y.coefs)));流(“maxdif = % g \ n”maxdif)如果(maxdif <公差),打破，结束%现在重申y = z;结束
Maxdif = 0.0184488 Maxdif = 0.000120467 Maxdif = 4.78115e-09
甘氨胆酸ax =;h = ax.Children;图例(h([10 9 5 4])，．．．｛= 1的初始猜想(x^2-1)“迭代”．．．sprintf ('对= %.3f的初步猜测'ε)．．．“迭代”}，“位置”，“西北”）;
非常に小さいイプシロン
εが非常に小さい場合は,毎回εを3で除算してこれらの計算を繰り返すだけです。
为e = 1:4节= newknt(z, ninterv+1);休息= knt2brk(节);节= augknt(休息,4,2);n =长度(节)- k;ninterv =长度(减免)1;temp =((优惠(2:ninterv + 1) +优惠(1:ninterv)) / 2);Temp = Temp ([1 1]，:) + gauss*diff(breaks); / /去掉所有的空colsites =临时(:)。';colsites网站= [0,1];Colmat = spcol(结点，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”) = /3;而1 vtau = fnval(y,colsites);Weights = [0 1 0;(2 * vtau。0 (n - 1) repmat(ε- 1));1 0 0];colloc = 0 (n, n);为j = 1: n colloc (j:) =重量(j:) * colmat (3 * (j - 1) + (1:3):);结束Coefs = colloc\[0 vtau.]* vtau + 1 0]。';z = spmak(节,系数。');fnplt (z,“b”）;maxdif = max (max (abs (z.coefs-y.coefs)));如果(maxdif <公差),打破，结束%现在重申y = z;结束结束甘氨胆酸ax =;h = ax.Children;图例(h([30 29 25 24])，．．．｛= 1的初始猜想(x^2-1)“迭代”．．．' = .1/3^j, j=1:5的初始猜测'．．．“迭代”}，“位置”，“西北”）;
最も小さいイプシロンに使用したブレークのプロット
次に,休息时间の最終的な分布を示します。この場合はnewkntがうまく機能していることがわかります。
休息= fnbrk (fn2fm (z,“页”)，“b”）;bb = repmat(休息,3,1);cc = repmat([0, 1,南],1,长度(休息);情节(bb (:), cc (:),“r”）;持有从甘氨胆酸ax =;h = ax.Children;Legend (h([31 30 26 25 1])，．．．｛= 1的初始猜想(x^2-1)“迭代”．．．' = .1/3^j, j=1:5的初始猜测'．．．“迭代”' = .1/3^5的断点'}，“位置”，“西北”）;
最も小さいイプシロンの残差のプロット
次の颂歌を解こうとしていることを思い出してください。
在[0..1]上D^2g(x) + (g(x))^2 = 1
チェックとして,最も小さいイプシロンについて計算した解の残差を計算およびプロットします。これも満足のいく結果に見えます。
xx = linspace (0, 1201);*fnval(fnder(z,2)，xx) - (fnval(z,xx)).^2)最小的计算解的残差）;
你也可以从以下列表中选择一个网站:
