理解Matlab中的DEL2函数，以便在C ++中进行编码 [英] Understanding DEL2 function in Matlab in order to code it in C++

问题描述

为了在c ++中编写DEL2 matlab函数，我需要理解算法。我设法编写了矩阵中不在边界或边缘的元素的函数。
我已经看到了几个关于它的主题，并通过键入编辑del2或键入del2来阅读MATLAB代码，但我不明白为获得边界和边缘所做的计算。

任何帮助将不胜感激，谢谢。 你想近似于你只知道你在一个点的右边（或左边）的价值。
为了有一个二阶近似，需要3个方程（基本的泰勒展开）：
$ b $ u（i + 1）= u（i）+ （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） （9/2）h ^ 2 u''+（27/6）h ^ 3（b + b）b
$ b u（i + 3）= u u'''+ O（h ^ 4）

解答u''给出（1）

：$ b $ （i + 1）+ 4（i + 2） - u（i + 3）+ 2 u（i）+ O（h ^ 4） ）

为了得到拉普拉斯，你需要用这个边界替换传统的公式。

（i = 0，j）= [-5 u（i + 1，j））的例子，其中i = 0 （i，j + 1）+ u（i，j-1）-2u（i，j + 1）+ u（i + 2，j） j）] / h ^ 2

编辑澄清：

x和y方向上的二阶导数之和。你可以用下面的公式计算二阶导数：（2）（）（）（）（）（）（） 1） - 2u（i））/ h ^ 2

如果你有u（i + 1）和u（i-1）。如果i = 0或i = imax，则可以使用我写的第一个公式计算导数（注意，由于二阶导数的扼制，如果i = imax，您可以用ik替换i + k）） 。同样适用于y（j）方向：

在边上，可以混合公式（1）和（ 2）：

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

在角落中，您只需在两个方向上使用（1）两次即可。

del2（u （i，j + 2） - u（i，j + 3）+ 2 u（i，j）+ （i，j + 1）+ 4 u（i + 2，j）-u（i + 3，j）+ 2 u（i，j）] / h ^ 2

Del2是二阶离散拉普拉斯算子，也就是说，它允许逼近连续函数的拉普拉斯算子，给定其在正方形笛卡尔网格上的值NxN，其中两个相邻n odes是 h 。

h ^ 2只是一个恒定的尺寸因子，您可以通过设置h ^ 2 = 4。

例如，如果要计算（0，L）上的u（x，y）的真实拉普拉斯算子）x（0，L）square，你要做的是在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 ^ 2的值放在矩阵A中。然后你将有
ans = 4 * del2（A）/ h ^ 2，其中h = L /（N-1）。

del2将返回如果你的起始函数是线性的或二次的（x ^ 2 + y ^ 2罚款，x ^ 3 + y ^ 3罚款），那么连续的拉普拉斯算子的精确值。如果函数不是线性的或二次的，那么结果将会更准确，因为您使用的点越多（即在极限h - > 0）。

我希望这是更多清楚，请注意，我使用基于0的索引访问矩阵（C / C ++数组样式），而matlab使用基于1的数据。

in order to code the DEL2 matlab function in c++ I need to understand the algorithm. I've managed to code the function for elements of the matrix that are not on the borders or the edges. I've seen several topics about it and read the MATLAB code by typing "edit del2" or "type del2" but I don't understand the calculations that are made to obtain the borders and the edges.

Any help would be appreciated, thanks.

解决方案

You want to approximate u'' knowing only the value of u on the right (or the left) of a point. In order to have a second order approximation, you need 3 equations (basic taylor expansion):

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)

Solving for u'' gives (1):

h^2 u'' = -5 u(i+1) + 4 u(i+2) - u(i+3) + 2 u(i) +O(h^4)

To get the laplacian you need to replace the traditional formula with this one on the borders.

For example where "i = 0" you'll have:

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

EDIT clarifications:

The laplacian is the sum of the 2nd derivatives in the x and in the y directions. You can calculate the second derivative with the formula (2)

u'' = (u(i+1) + u(i-1) - 2u(i))/h^2

if you have both u(i+1) and u(i-1). If i=0 or i=imax you can use the first formula I wrote to compute the derivatives (notice that due to the simmetry of the 2nd derivative, if i = imax you can just replace "i+k" with "i-k"). The same applies for the y (j) direction:

On the edges you can mix up the formulas (1) and (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

And on the corners you'll just use (1) two times for both directions.

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 is the 2nd order discrete laplacian, i.e. it permits to approximate the laplacian of a real continuous function given its values on a square cartesian grid NxN where the distance between two adjacent nodes is h.

h^2 is just a constant dimensional-factor, you can get the matlab implementation from these formulas by setting h^2 = 4.

For example, if you want to compute the real laplacian of u(x,y) on the (0,L) x (0,L) square, what you do is writing down the values of this function on an NxN cartesian grid, i.e. you calculate 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)) etc. and you put down these N^2 values in a matrix A.

Then you'll have ans = 4*del2(A)/h^2, where h = L/(N-1).

del2 will return the exact value of the continuous laplacian if your starting function is linear or quadratic (x^2+y^2 fine, x^3 + y^3 not fine). If the function is not linear nor quadratic, the result will be more accurate the more points you use (i.e. in the limit h -> 0)

I hope this is more clear, notice that i used 0-based indices for accessing matrix (C/C++ array style), while matlab uses 1-based.

这篇关于理解Matlab中的DEL2函数，以便在C ++中进行编码的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！