2

C++でDEL2matlab関数をコーディングするには、アルゴリズムを理解する必要があります。境界やエッジにない行列の要素の関数をコーディングすることができました。それに関するいくつかのトピックを見て、「editdel2」または「typedel2」と入力してMATLABコードを読みましたが、境界線とエッジを取得するために行われる計算がわかりません。

助けていただければ幸いです、ありがとう。

4

4 に答える 4

5

ポイントの右側（または左側）のuの値だけを知って、u''を近似したいとします。2次近似を行うには、次の3つの方程式が必要です（基本的なテイラー展開）。

u（i + 1）= u（i）+ h u'+（1/2）h ^ 2 u''+（1/6）h ^ 3 u''' + O（h ^ 4）

u（i + 2）= u（i）+ 2 h u'+（4/2）h ^ 2 u''+（8/6）h ^ 3 u''' + O（h ^ 4）

u（i + 3）= u（i）+ 3 h u'+（9/2）h ^ 2 u''+（27/6）h ^ 3 u''' + O（h ^ 4）

u''を解くと（1）が得られます：

h ^ 2 u'' = -5 u（i + 1）+ 4 u（i + 2）-u（i + 3）+ 2 u（i）+ O（h ^ 4）

ラプラシアンを取得するには、境界で従来の式をこれに置き換える必要があります。

たとえば、「i = 0」の場合、次のようになります。

del2（u）（i = 0、j）= [-5 u（i + 1、j）+ 4 u（i + 2、j）-u（i + 3、j）+ 2 u（i、j） + u（i、j + 1）+ u（i、j-1）-2u（i、j）] / h ^ 2

説明を編集する：

ラプラシアンは、x方向とy方向の2次導関数の合計です。式（2）を使用して2階導関数を計算できます。

u'' =（u（i + 1）+ u（i-1）-2u（i））/ h ^ 2

u（i + 1）とu（i-1）の両方がある場合。i=0またはi=imaxの場合、導関数を計算するために私が書いた最初の式を使用できます（2次導関数のシミュレーションにより、i = imaxの場合、「i +k」を「ik」に置き換えることができます） 。同じことがy（j）方向にも当てはまります。

エッジでは、式（1）と（2）を混同することができます。

del2（u）（i = imax、j）= [-5 u（i-1、j）+ 4 u（i-2、j）-u（i-3、j）+ 2 u（i、j） + u（i、j + 1）+ u（i、j-1）-2u（i、j）] / h ^ 2

del2（u）（i、j = 0）= [-5 u（i、j + 1）+ 4 u（i、j + 2）-u（i、j + 3）+ 2 u（i、j） + u（i + 1、j）+ u（i-1、j）-2u（i、j）] / h ^ 2

del2（u）（i、j = jmax）= [-5 u（i、j-1）+ 4 u（i、j-2）-u（i、j-3）+ 2 u（i、j） + u（i + 1、j）+ u（i-1、j）-2u（i、j）] / h ^ 2

そして、コーナーでは、（1）を両方向に2回使用します。

del2（u）（i = 0、j = 0）= [-5 u（i、j + 1）+ 4 u（i、j + 2）-u（i、j + 3）+ 2 u（i、 j）+ -5 u（i、j + 1）+ 4 u（i + 2、j）-u（i + 3、j）+ 2 u（i、j）] / h ^ 2

Del2は、2次の離散ラプラシアンです。つまり、2つの隣接するノード間の距離がhである正方形のデカルトグリッドNxNでの値が与えられると、実際の連続関数のラプラシアンを近似できます。

h ^ 2は単なる一定の次元係数であり、h ^ 2 = 4を設定することにより、これらの式からmatlabの実装を取得できます。

たとえば、（0、L）x（0、L）の正方形でu（x、y）の実際のラプラシアンを計算する場合は、この関数の値をNxNデカルトグリッドに書き留めます。つまり、u（0,0）、u（L /（N-1）、0）、u（2L /（N-1）、0）... u（（N-1）L /（N- 1）= L、0）... u（0、L /（N-1））、u（L /（N-1）、L /（N-1））などとあなたはこれらのN^を置きます行列Aの2つの値。

次に、ans = 4 * del2（A）/ h ^ 2になります。ここで、h = L /（N-1）です。

開始関数が線形または2次の場合、del2は連続ラプラシアンの正確な値を返します（x ^ 2 + y ^ 2は細かい、x ^ 3 + y ^ 3は正しくない）。関数が線形でも二次でもない場合、使用するポイントが多いほど結果はより正確になります（つまり、限界h-> 0）

これがより明確になることを願っています。matlabは1ベースを使用しているのに対し、マトリックス（C / C ++配列スタイル）にアクセスするために0ベースのインデックスを使用していることに注意してください。

于 2012-11-12T11:36:33.540 に答える
1

MatLabのDEL2は離散ラプラス演算子を表します。ここでいくつかの情報を見つけることができます。

エッジに関する主な点は、マトリックスの内部の要素には4つの隣接要素があり、エッジとコーナーの要素にはそれぞれ3つまたは2つの隣接要素があることです。したがって、コーナーとエッジを同じ方法で計算しますが、使用する要素は少なくなります。

于 2012-11-12T11:35:47.940 に答える
0

これは、上記のアイデアを実装するMATLABの「del2（）」演算子を複製するFortran90で作成したモジュールです。少なくとも4x4以上のアレイでのみ機能します。実行するとうまくいくので、他の人が自分で時間を無駄にしないように投稿したいと思いました。

module del2_mod
implicit none
real, private                       :: pi
integer, private                    :: nr, nc, i, j, k
contains
! nr is number of rows in array, while nc is the number of columns in the array.
!!---------------------------------------------------------- 

subroutine del2(in, out)
real, dimension(:,:)            :: in, out
real, dimension(nr,nc)          :: interior, left, right, top, bottom, ul_corner, br_corner, disp
integer                         :: i, j
real                            :: h, ul, ur, bl, br
! Zero out internal arrays
out = 0.0; interior=0.0; left = 0.0;  right = 0.0;  top = 0.0;  bottom = 0.0;  ul_corner = 0.0; br_corner = 0.0;
h=2.0

! Interior Points
do j=1,nc
    do i=1,nr
    ! Interior Point Calculations
    if( j>1 .and. j<nc .and. i>1 .and. i<nr )then
        interior(i,j) = ((in(i-1,j) + in(i+1,j) + in(i,j-1) + in(i,j+1)) - 4*in(i,j) )/(h**2)
    end if
    ! Boundary Conditions for Left and Right edges
    left(i,1) = (-5.0*in(i,2) + 4.0*in(i,3) - in(i,4) + 2.0*in(i,1) + in(i+1,1) + in(i-1,1) - 2.0*in(i,1) )/(h**2)
    right(i,nc) = (-5.0*in(i,nc-1) + 4.0*in(i,nc-2) - in(i,nc-3) + 2.0*in(i,nc) + in(i+1,nc) + in(i-1,nc) - 2.0*in(i,nc) )/(h**2)
    end do
! Boundary Conditions for Top and Bottom edges
top(1,j) = (-5.0*in(2,j) + 4.0*in(3,j) - in(4,j) + 2.0*in(1,j) + in(1,j+1) + in(1,j-1) - 2.0*in(1,j) )/(h**2)
bottom(nr,j) = (-5.0*in(nr-1,j) + 4.0*in(nr-2,j) - in(nr-3,j) + 2.0*in(nr,j) + in(nr,j+1) + in(nr,j-1) - 2.0*in(nr,j) )/(h**2)
end do
out = interior + left + right + top + bottom 
! Calculate BC for the corners
ul = (-5.0*in(1,2) + 4.0*in(1,3) - in(1,4) + 2.0*in(1,1) - 5.0*in(2,1) + 4.0*in(3,1) - in(4,1) + 2.0*in(1,1))/(h**2)
br = (-5.0*in(nr,nc-1) + 4.0*in(nr,nc-2) - in(nr,nc-3) + 2.0*in(nr,nc) - 5.0*in(nr-1,nc) + 4.0*in(nr-2,nc) - in(nr-3,nc) + 2.0*in(nr,nc))/(h**2)
bl = (-5.0*in(nr,2) + 4.0*in(nr,3) - in(nr,4) + 2.0*in(nr,1) - 5.0*in(nr-1,1) + 4.0*in(nr-2,1) - in(nr-3,1) + 2.0*in(nr,1))/(h**2)
ur = (-5.0*in(1,nc-1) + 4.0*in(1,nc-2) - in(1,nc-3) + 2.0*in(1,nc) - 5.0*in(2,nc) + 4.0*in(3,nc) - in(4,nc) + 2.0*in(1,nc))/(h**2)
! Apply BC for the corners
out(1,1)=ul
out(1,nc)=ur
out(nr,1)=bl
out(nr,nc)=br
end subroutine

end module
于 2013-08-02T00:39:29.883 に答える
0

難しいです！私はそれを理解してJavaで実装するのに数時間を無駄にしました。

ここにあります：https ：//gist.github.com/emersonmoretto/dec8f7125c032775da0d

テストされ、元の関数DEL2（Matlab）と比較されました

sbabbiの応答にタイプミスが見つかりました：

del2(u) (i=0,j=0) = [-5 u(i,j+1) + 4 u(i,j+2) - u(i,j+3) + 2 u(i,j) + -5 u(i,j+1) + 4 u(i+2,j) - u(i+3,j) + 2 u(i,j)]/h^2

は

del2(u) (i=0,j=0) = [-5 u(i,j+1) + 4 u(i,j+2) - u(i,j+3) + 2 u(i,j) + -5 u(i+1,j) + 4 u(i+2,j) - u(i+3,j) + 2 u(i,j)]/h^2
于 2014-09-03T04:34:22.723 に答える