在 MATLAB 中使用移位矩阵求数值导数

Question

我们被分配使用给定步骤找到数值一阶和二阶导数：

1)定义两个数组，x和y=f(x)(使用任意函数)

2) 将微分算子定义为矩阵（例如，对于 f' ，定义矩阵 D= (1/2h)(U-L) 其中 U 是 Uppershift 矩阵，L 是 Lowershift 矩阵。

3)微分算子矩阵乘以函数数组求导

我写了这段代码：

%QUESTION 3- DIFFERENTIAL OPERATOR
h=2;
x = 2:h:8
y = x.^2 %the chosen function


n=length(x);
shift = 1; 
U = diag(ones(n-abs(shift),1),shift);
s = -1;
L= diag(ones(n-abs(shift),1),s);
% % the code above creates the upper and lower shift matrix
%
D= ((U-L))/(2*h) %differential operator 
d = y*D %approximates the first derivative of each element in vector -x
d2 = d*D %approximates the second derivative of each element in vector -x

为了创建移位矩阵我使用了

现在我得到了这个作为解决方案：

>> mat2q3

x =

     2     4     6     8


y =

     4    16    36    64


D =

         0    0.2500         0         0
   -0.2500         0    0.2500         0
         0   -0.2500         0    0.2500
         0         0   -0.2500         0


d =

    -4    -8   -12     9


d2 =

    2.0000    2.0000   -4.2500   -3.0000

请告诉我我哪里做错了，或者我该如何改进。

Answer 1

这段代码的大部分功能都在做它打算做的事情：计算导数的数值近似值。但是，您计算导数的方式有点不正确。基本上，对于数组中的每个点，您要用右边的点减去左边的点，然后除以 2*h。如果你想这样做，你需要 post- 乘以向量 y 并将其转换为列向量。

转置 y 然后取 D 并乘以 y:

>> d = D*y.'

d =

     4
     8
    12
    -9

但是，我想指出第一个和最后一个条目没有任何意义，因为您的 D 矩阵：

D =

         0    0.2500         0         0
   -0.2500         0    0.2500         0
         0   -0.2500         0    0.2500
         0         0   -0.2500         0

您在这里计算的是数值导数的中心差，它需要在您计算的点的左侧和右侧有一个点大批。基本上，对于第一点，您只有右侧的点，而当您越界时，左侧的点是未定义的。同样对于最后一点，你只有左边的一个点，右边的点也不存在，因为它不存在。

此处的数值导数矩阵假设超出范围的点为0。事实证明你在这一点上获得了正确的导数，这是侥幸（即2*x在x = 1 是 2)。原因是因为 x = 0 处的点确实使 x.^2 的输出等于 0，所以在数组中省略这个点就好像它在你的点数组中一样开始。

二阶导数同样需要：

>> d2 = D*d

d2 =

    2.0000
    2.0000
   -4.2500
   -3.0000

但是，请记住，对于您之前的一阶导数结果 d，最后一项是垃圾，因此如果您计算该结果的导数，您将使用不正确的最右项，当您查看 d 的第三个条目，因此 d2 的最后两个条目和前两个条目都是垃圾。恰好 d2 的第二个条目是正确的。这归因于 x = 0 实际上并不存在于先前的结果中，但是假设您的点数组中的值越界为 0，我们确实得到了正确的结果。

您应该在更长的序列上尝试此操作。例如，尝试这样做：

h = 2;
x = 2:h:20;
y = x.^2;

这定义了一个从 2 到 20 的序列，步长为 h = 2。这就是 x 和 d:

的结果

x =

    2     4     6     8    10    12    14    16    18    20

d =

     4
     8
    12
    16
    20
    24
    28
    32
    36
   -81

如您所见，第一个和最后一个条目没有意义，因为我们没有第一个条目的最左边点和第二个条目的最右边点来考虑。然而，如前所述，侥幸的第一点给了我们正确的结果。一般来说，第一个和最后一个条目不应该给你正确的结果，因为我们没有左点或右点来解释导数。对于其余的点，您可以看到我们正在正确计算导数，即 2*x。

我们来看看二阶导数，d2:

碰巧前两个元素是正确的，即使它们在理论上不应该是正确的，但我们侥幸得到了它。但是，最后两个元素不正确。

你需要从中得到的是记住一件事。请记住，当使用中心差进行数值导数时，您必须 trim 前 n 个元素和最后 n 个元素，其中 n 是导数的顺序你在考虑。

在 MATLAB 中使用移位矩阵求数值导数

Finding Numerical Derivative using a shift matrix in MATLAB

matlab

matrix

numerical-methods