matlab - Streamslice 函数可能存在错误，或者只是在 MATLAB 中误用了它？

Question

问题很简单也很直接:我想绘制条形磁铁(等于螺线管磁场)的磁流，计算公式为 these equations 。我的 MATLAB 版本是 R2013a。

完成这项工作的函数，Magnetic_field(它返回一个具有坐标(r，z)和场(Br，Bz)的结构，B(r ,z) = Br(r,z)r + Bz(r,z)z 已计算)，在本题末尾。

然后我打电话

>> t = magnetic_field(2e-2, 16e-2, 20, -6e-2, 6e-2, 0.5e-2, -18.2e-2, 18.2e-2, 1e-2);
>> streamslice(t.r, t.z, t.Br, t.Bz)

请注意，t.r 和 t.z 是一个 meshgrid 结果，其中 Br 和 Bz 进行计算(请参阅函数 Magnetic_field 中的 for 循环)。因此，所有这四个矩阵都具有相同的大小。

MATLAB 给出错误:

Attempted to access startgrid(NaN,13); index must be a positive integer or logical.

Error in streamslice>nicestreams (line 324) startgrid(rr,cc)=1;

Error in streamslice (line 134) [vertsout, avertsout] = nicestreams(x,y,u,v,density,arrows);

奇怪的是:如果我更改网格，streamslice 不会给出错误:

>> t = magnetic_field(2e-2, 16e-2, 20, -2, 2, 0.1, -2, 2, 0.1);

尽管它产生了一个非常奇怪的结果(因为我认为该解决方案仅针对正 r 定义)。

我去检查了 MATLAB 示例数据，例如 wind 数据:

>> w = load('wind');
>> streamslice(w.x(:,:,5), w.y(:,:,5), w.u(:,:,5), w.v(:,:,5))

它绘制正确。但是，如果我从 5 更改为 2、3、10、12 > 或 15，MATLAB 显示与我遇到的错误非常相似的错误:

Attempted to access startgrid(2,0); index must be a positive integer or logical.

Error in streamslice>nicestreams (line 324) startgrid(rr,cc)=1;

Error in streamslice (line 134) [vertsout, avertsout] = nicestreams(x,y,u,v,density,arrows);

不同之处在于，对于示例数据，错误在于 streamslice 尝试访问 0 索引，而在我的例子中，它尝试访问 >NaN 索引。

因此，我保留两种可能性:

1) 我的数据(即变量 Br 和 Bz)不具有 streamslice 函数所需的结构 2)这是streamslice函数中的一个错误

你觉得怎么样？

有什么建议吗？

function res = magnetic_field(a, L, I, rMin, rMax, dr, zMin, zMax, dz)    
    Nr = (rMax - rMin) ./ dr + 1;
    Nz = (zMax - zMin) ./ dz + 1;

    % calculating dependent constants
    mu0div2pi = 2e-7;
    Lm1 = 1 / L; % in 1/m
    Ldiv2 = L / 2;
    constZ = Lm1 .* I * mu0div2pi / 4;
    constR = Lm1 .* I * mu0div2pi;

    % defining functions
    kSqr = @(ar4, ar2, zeta) ar4 ./ (ar2 + zeta .^ 2);
    hSqr = @(ar4, ar2) ar4 ./ ar2;
    zetaPlus = @(z, Ldiv2) z + Ldiv2;
    zetaMinus = @(z, Ldiv2) z - Ldiv2;
    IntBr = @(kSqr, eK) ((kSqr - 2) .* eK + 2 .* ellipticE(kSqr)) ./ sqrt(kSqr);
    IntBz = @(kSqr, hSqr, r, zeta, a, eK) zeta .* sqrt(kSqr) .* (eK + ellipticPi(hSqr, kSqr) .* (a - r) ./ (a + r));

    % creating grid
    [r, z] = meshgrid(linspace(rMin, rMax, Nr), linspace(zMin, zMax, Nz));
    [m, n] = size(r);
    NN = m * n;

    Br = zeros(m,n);
    Bz = zeros(m,n);
    for i = 1:NN
        ar4 = 4 .* a .* r(i);
        ar2 = (a + r(i)) .^ 2;
        hh = hSqr(ar4, ar2);
        zP = zetaPlus(z(i), Ldiv2);
        zM = zetaMinus(z(i), Ldiv2);
        kkPlus = kSqr(ar4, ar2, zP);
        kkMinus = kSqr(ar4, ar2, zM);
        eKPlus = ellipticK(kkPlus);
        eKMinus = ellipticK(kkMinus);
        Br(i) = constR .* sqrt(a ./ r(i)) .* (IntBr(kkPlus, eKPlus) - IntBr(kkMinus, eKMinus));
        Bz(i) = constZ .* ( IntBz(kkPlus, hh, r(i), zP, a, eKPlus) - IntBz(kkMinus, hh, r(i), zM, a, eKMinus) ) ./ sqrt(a .* r(i));
    end
    res.a = a;
    res.L = L;
    res.I = I;
    res.r = r;
    res.z = z;
    res.Br = Br;
    res.Bz = Bz;
end

Answer 1

这似乎确实是 streamslice 中的错误正如 @user664303 提到的，你应该报告它。

---

w = load('wind');
i = 2;
streamslice(w.x(:,:,i), w.y(:,:,i), w.u(:,:,i), w.v(:,:,i))

使用wind按照您的指示(我使用的是 MATLAB R2013b)数据集，该错误似乎是由舍入误差引起的；值xc-xmin意外地返回一个接近机器精度的非常小的负数，该负数又变成 -1调用floor功能。该结果(加一)稍后用于索引，从而导致错误。

好消息是我相信这里有一个简单的修复方法；将调用替换为 floor调用fix相反(没有双关语的意思)。这应该执行向零舍入，而不是向负无穷大舍入，这在正数的情况下(在这种情况下似乎是预期的)是等效的，同时捕获轻微的舍入错误，当应该出现零时，可能会导致非常小的负数，并进行纠正。

我必须在 streamslice.m 内的六个不同位置执行此更改。文件(只需将所有出现的 floor 替换为 fix )。这是更改的补丁文件:

--- old/streamslice.m
+++ new/streamslice.m
@@ -321,8 +321,8 @@
                     yc = vv(j,2);
-                    cc = floor( (xc-xmin)*ixrangecs )+1;
-                    rr = floor( (yc-ymin)*iyrangers )+1;
+                    cc = fix( (xc-xmin)*ixrangecs )+1;
+                    rr = fix( (yc-ymin)*iyrangers )+1;
                     startgrid(rr,cc)=1;

-                    cc = floor( (xc-xmin)*ixrangece )+1;
-                    rr = floor( (yc-ymin)*iyrangere )+1;
+                    cc = fix( (xc-xmin)*ixrangece )+1;
+                    rr = fix( (yc-ymin)*iyrangere )+1;
                     if endgrid(rr,cc)==1
@@ -336,4 +336,4 @@
                     if arrows
-                        cc = floor( (xc-xmin)*ixrangeca )+1;
-                        rr = floor( (yc-ymin)*iyrangera )+1;
+                        cc = fix( (xc-xmin)*ixrangeca )+1;
+                        rr = fix( (yc-ymin)*iyrangera )+1;
                         if j>1 && arrowgrid(rr,cc)==0

比较R2013a和R2013b的功能，文件是相同的。所以我认为这个解决方法也应该适用于该版本。请在实际数据上尝试一下，看看是否有效。

---

编辑:

正如您所发现的，streamslice 中似乎存在第二个错误。功能。它仅在数据包含 NaN 值时影响代码。据我所知，解决方法如下:

--- old/streamslice.m
+++ new/streamslice.m
@@ -302,7 +302,7 @@
             for q = 1:2
                 vv = verts{q};

-                nanpos = find(isnan(vv));
+                nanpos = find(any(isnan(vv),2));
                 if ~isempty(nanpos)
                     if nanpos(1)==1
                         vv = [];

值得注意的是，使用您的函数生成数据:

t = magnetic_field(2e-2, 16e-2, 20, -6e-2, 6e-2, 0.5e-2, -18.2e-2, 18.2e-2, 1e-2);

在 R2013a 和 R2013b 之间给出不同的结果(特别是使用 R2013a，t.Bz 有一列包含 R2013b 生成的数据中不存在的所有 NaN/Inf 值)...这解释了为什么 Nan-bug 没有出现在 R2013b 中。