MATLAB 使用 ode45 对波函数进行动画处理会引发容差错误答案

【问题标题】：MATLAB using ode45 to animate wave function throws tolerance errorMATLAB 使用 ode45 对波函数进行动画处理会引发容差错误
【发布时间】：2020-04-12 23:08:22
【问题描述】：

我正在制作 MATLAB 中的波函数动画。到目前为止，我还没有接触到动画的一部分，因为我的 ode45 引发了错误。我为此使用 de Korteweg-de Vries 方程（可以在这里找到：https://en.wikipedia.org/wiki/Korteweg%E2%80%93de_Vries_equation）。

现在我的代码如下：

function[t, N, u] = wave(u0, L, T, S, N)

%First we create a vector t which stores the time entries, a vector x which
%stores the location entries and an h which is the location step size for
%the Korteweg-De Vries function.
time = linspace(0,T,S);
x = linspace(-L,L,N);
h = x(2)-x(1);

U0=u0(x);

options=odeset('RelTol',1e-13,'AbsTol',1e-13);
[t,u] = ode45(@kdv,time,U0,options);

    function dudt = kdv(t,u)

    d2udx2 = ([u(N)-2*u(1)+u(2); diff(u,2); u(N-1)-2*u(N)+u(1)] ./ h^2);
    total = d2udx2 + 3.*u.^2;

    diff_total = diff(total);
    diff_total = [diff_total(end); diff_total(1:end)];

    dudt = -[diff_total(2)-diff_total(N-1); diff(diff_total,2); diff_total(N-1)+diff_total(2)] ./ (2*h);

    end
end

现在，当我设置 f=@(x)(2/2).*(1./cosh(sqrt(2).*x./2).^2) 然后使用 [t,N,u]=wave(f, 20, 10, 200, 200); 调用函数时，我收到以下错误：

警告：在 t=6.520003e-02 失败。如果不将步长减小到低于在时间 t 允许的最小值 (2.220446e-16)。在 ode45 中（第 360 行）在波浪中（第 23 行）

返回的t 和u 的大小分别为16x1 和16x200，但我认为如果没有出现警告，它们会扩展为200x1 和200x200。

任何有关如何解决此问题的提示将不胜感激。

【问题讨论】：

使用diff_total将迭代差异增加到订单5的原因是什么？为什么最终的导数向量存在二阶差分？为什么不通过h 的适当幂除将迭代的差异变成适当的差商？
Korteweg-de vries 方程表明 du/dt = -d/dx * (d^2u/dx^2 + 3 u^2)。因此，我认为我需要二阶差异，才能将那部分放在括号之间，对吧？
忘记我的最后一句话，正确的划分就在那里。 d2udx2 和 total 构造正确。移动其他评论以回答。

标签： matlab warnings integration ode

【解决方案1】：

你做的差异太大，构造diff_total的中间步骤是多余的。

d2udx2 和 total 构造正确。您可以将第一个构造缩短为
```
d2udx2 = diff([u(N);  u;  u(1) ], 2)./h^2`;
```
下一步构造 diff_total 是多余的，与 dudt 中的差异构造不同。 diff_total 构造缺少 h 的除法，由于某种未知原因，您在 dudt 中间有差异顺序 2。
外部x 微分最好通过卷积来获得中心差商，就像在外部元素中一样，
```
dudt = -conv([ total(N); total; total(1)], [1; 0; -1], 'valid')/(2*h)
```

所有这些都给出了导数过程

function dudt = kdv(t,u,h)
    d2udx2 = diff([u(end);  u;  u(1) ], 2)./h^2;
    total = d2udx2 + 3*u.^2;
    dudt = -conv( [ total(end); total; total(1)], [1; 0; -1], 'valid')./(2*h);
end

which (with N=80) 产生解决方案

正如预期的那样，它是以c=2 的速度传播的孤子波。

【讨论】：