1

関数 @f(x,y) が与えられ、MATLAB で特定の凸多角形に対するこの関数の積分を評価したいと考えています。多角形は必ずしも四角形であるとは限らないため、MATLAB の関数「dblquad」を使用できません。私が持っているポリゴンは、ベクトル X と Y で表される頂点のセットによって与えられます。つまり、頂点は (X(1),Y(1)),....,(X(n),Y(n) ）。使用できる関数またはメソッドはありますか?

4

1 に答える 1

2

秘訣は、ツールを使用して関心領域内に統合することです。三角領域で統合するためのいくつかのツールを作成しました。

% Define a function to integrate.
% This function takes an nx2 array, where each row
% contains a single point to evaluate the kernel at.
% This computes x^2 + y^2 at each point.
fun = @(xy) sum(xy.^2,2);

% define the domain as a triangulated polygon
% this tool uses ear clipping to do so.
sc = poly2tri([1 4 3 1],[1 3 5 4]);

% Gauss-Legendre integration over the 2-d domain
[integ,fev]= quadgsc(fun,sc,2)
integ =
       113.166666666667
fev =
     8

% the triangulated polygon...
plotsc(sc,'facecolor','none','markerfacecolor','r')
axis equal
grid on

ここに画像の説明を入力

その多角形ドメイン上のマッピング z(x,y) として、関数自体を視覚化できます。範囲フィールドが提供されると、シンプリシアル複体は 2-d (x,y) ドメインからの 2-1 マッピングに変わります。

sc2 = refinesc(sc,'max',.5);
sc2.range = fun(sc2.domain);
plotsc(sc2,'markerfacecolor','r')
grid on
view(17,12)

ここに画像の説明を入力

これは、関心のあるドメイン上の単純な多項式関数であるため、デフォルトの低次ガウス積分で十分でした。使用されるスキームは、三角形上のテンソル積形式のガウス ルジャンドル スキームであり、真に最適ではありませんが実行可能です。ガウス求積法の問題は、適応性がないことです。点の有限集合に対する多項式による暗黙的な近似に基づいて、推定値を計算します。

上記の見積もりでは、その見積もりを計算するために 8 つの関数評価が使用されました。カーネルは低次多項式であるため、完全に機能するはずです。問題は、それが正しい解決策であるかどうかを知る必要があるということです。これはガウス求積法の問題です。収束するように見えるまで高次スキームで問題を解決する以外に、答えが正しいかどうかを知る簡単な方法はありません。

重心で三角形ごとに 1 つのポイントを使用すると、間違った答えが得られますが、高次の推定値はすべて一致します。

[integ,fev]= quadgsc(fun,sc,1)
integ =
       107.777777777778
fev =
     2

[integ,fev]= quadgsc(fun,sc,3)
integ =
       113.166666666667

fev =
    18

[integ,fev]= quadgsc(fun,sc,4)
integ =
       113.166666666667
fev =
    32

quadgsc を作成した後、MATLAB の他のクワッド ツールと同じように機能する適応ソルバーを試す必要がありました。これは、三角形分割の適応改良を行い、解が安定していない三角形を探します。問題は、これらのツールを満足のいくように書き終えたことがないことです。三角領域上の立方体問題に使用できるさまざまな方法があります。quadrsc は、低次の解を実行してからそれを改良し、リチャードソン外挿を使用して、結果を比較します。差が大きすぎる三角形については、収束するまでさらに改良します。

例えば、

[integ,fev]= quadrsc(fun,sc)
integ =
          113.166666666667

fev =
    16

したがって、これは機能します。この問題は、より複雑なカーネルで発生します。この問題は、改良をいつ停止するかを知ることになり、多くの関数評価を使い果たす前に停止する必要があります。私はそれを完全に満足のいくように機能させることができなかったので、これらのツールを投稿したことはありません. ダイレクトメールを送ってくださった方にはツールボックスをお送りできます。zip ファイルは約 2.4 MB です。いつの日か、それらのツールを完成させられる日が来ることを願っています...

于 2013-02-05T11:46:06.093 に答える