MATLAB 求解曲线不适合使用 quiver 和 dsolve 的斜率场答案

【问题标题】：MATLAB Solution curves not fitting slope field using quiver and dsolveMATLAB 求解曲线不适合使用 quiver 和 dsolve 的斜率场
【发布时间】：2018-06-15 05:04:49
【问题描述】：

给定一个一阶微分方程...

dy/dx = 0.01y(1 - y)(y - 1/2)

我正在尝试创建一个斜率场并显示一些特定的解曲线，但由于某种原因，这些曲线似乎不能完美地拟合向量。

这是输出：如果您查看两条中间曲线，它们的斜率似乎并不跟随向量。这是情节本身的产物还是我做错了什么？

我花了一些时间让斜率场很好地显示出来，也许解曲线不适合的原因是因为点之间的间距，但我不确定。

代码如下：

    [T Y] = meshgrid(-200:10:200, -2:0.1:2);
    dY = 0.01.*Y.*(1-Y).*(Y-(0.5));
    dT = ones(size(dY));
    L = sqrt(dT.^2 + dY.^2);
    figure
    quiver(T, Y, dT./(L), dY./L, .6, 'AutoScale', 'off',                         'ShowArrowHead', 'off');
    axis tight
    hold on
    grid on
    axis([-200 200 -2 2])

    syms x y
    de = 'Dy = 0.01*y*(1-y)*(y-(0.5))';

    y1 = dsolve(de, 'y(0) = -1', 'x');
    y1 = expand(y1);

    y2 = dsolve(de, 'y(0) = 1/4', 'x');
    y2 = expand(y2);

    y3 = dsolve(de, 'y(0) = 3/4', 'x');
    y3 = expand(y3);

    y4 = dsolve(de, 'y(0) = 2', 'x');
    y4 = expand(y4);

    for fxn = [y1 y2 y3 y4]
        fplot(fxn, [-200 200], 'LineWidth', 2)
        hold on
    end

【问题讨论】：

是的，它看起来像是剧情的神器。中间的斜率使得分段中的高度差小于 1 个像素，解曲线也具有每 20 个分段 10-20 个像素的斜率。您可以在抗锯齿的污迹中看到一些正确的斜率。如果您要生成更高分辨率的图片，您应该会更好地看到斜坡段并不是真正水平的。

标签： matlab differential-equations graphing

【解决方案1】：

Quiver 在矩形网格 (https://www.mathworks.com/help/matlab/ref/quiver.html) 上绘制有明确定义的原点。而不是说大量带有不均匀虚线图的曲线，其中虚线将匹配该点的值。

在中心，线条几乎是水平的。因此，看起来好像你的紫色和橙色曲线的低斜率与颤动相交并且不遵循假设的线。如果您在曲线附近放大，您会看到差异，并且两条曲线实际上都正确地遵循斜率。

很遗憾，您无法通过使用 quiver 的更密集的网格来解决此问题。假设您将 10 倍的点添加到 T 和 Y - 所有的图形都将是一片蓝色的海洋，但在紫色线附近放大仍然会显示同样的问题。

现在，您可以在中心区域手动绘制“类似颤动”的图，例如在 Y 上从 0 到 1。选择起始 Y 和 X 点，可能与您现在的相同。画出与相同箭袋相对应的第一条线。现在，不要让其他向量的后续 Y 点保持相同，而是计算如果直线继续下去它们必须在的位置 - Ynext = Y+dY*(T(i)-T(i-1))。为下一轮点的修改后的 Y 点计算相同的 dY、dT 和 L，然后重复直到到达边界的末端。

说出这段代码（我仍然需要修正长度和校正偏移，但总体思路如下）：

Ycentral = 0:0.1:1;
Ynext = Ycentral;
numX = length(x);
for i = 1 : numX   
    [T2, Y2] = meshgrid(x(i), Ynext);
    dY = 0.01 * Ynext .* (1-Ynext).*(Ynext - 0.5);
    dT = ones(size(dY));
    L = sqrt(dT.^2 + dY.^2);    

    quiver(T2', Y2', dT./(L), dY./L, .6, 'AutoScale', 'off', 'ShowArrowHead', 'off');

    Ynext = Ynext + dY;
end

【讨论】：