MATLAB中ODE解导数的错误值

Question

用这里 (见，ode4)描述的经典Runge方法求解非线性ODE的柯西问题。在运行这部分之后

q = 10;
T = 5;
N = @(g1, g2) g1^2 + sin(g2);

t = linspace(0, T, q);
GCC = nonlinearGreenCC(N, 1, T);
di = diff(GCC);

我评估GCC(1)和di(1)。而GCC(1) = 0如预期，di(1) = 1.6e-05。不明白为什么，因为一阶导数上的柯西条件是1。如何修正不准确/错误？任何帮助都是非常感谢的。

函数nonlinearGreenCC如下：

function G = nonlinearGreenCC(N, a0, T)

h = .0001;
eqGreen = @(t, g)[g(2); - N(g(1), g(2))];
Cc = [0, a0];
G = ode4(eqGreen, 0, h, T, Cc);

end

Answer 1

GCC是一对值的数组，从更精确地与以下内容进行集成开始

[  0.00000000e+00,   1.00000000e+00],
[  9.99957927e-05,   9.99915855e-01],
[  1.99983171e-04,   9.99831715e-01],
[  2.99962136e-04,   9.99747579e-01],
[  3.99932687e-04,   9.99663448e-01],

对于确切的值，但是使用ode4函数的结果开始时如下所示

 0                         1
 1.66652642851402e-05      1.00001666526429
-1.40231141152723e-05      0.999985976885885
 1.66638620438405e-05      1.00001666386204
-1.40217119017434e-05      0.999985978288098

它不包含时间间隔的抽样，实际上，在您称为diff的级别上，h的值是未知的。diff不可能计算出您所期望的差值商。如果您阅读了文档，您应该会发现它只计算了差异。第一个差异只返回第一个值1.66652642851402e-05。

导致ode4算法产生奇怪结果的直接错误是，eqGreen生成应该返回行向量的列向量。由于初始值是行向量，在ode4中添加行向量和列向量的结果会产生一个2x2矩阵，该矩阵作为2行添加到结果中，结果会产生令人困惑的结果。两行向量都将结果以交替值和导数的形式放入一行。变到

eqGreen = @(t, g)[g(2), - N(g(1), g(2))];
Cc = [0, a0];

结果是

GCC(1:5,:) =

                     0                     1
  9.99957927208439e-05     0.999915855174488
  0.000199983171186392     0.999831714893813
  0.000299962135851046     0.999747579156327
  0.000399932687169042     0.999663447960381

di(1:5,:)=

  9.99957927208439e-05  -8.41448255122224e-05
  9.99873784655483e-05  -8.41402806746050e-05
  9.99789646646540e-05  -8.41357374861129e-05
  9.99705513179961e-05  -8.41311959463020e-05
  9.99621384254097e-05  -8.41266560547282e-05

如果用h=1e-4除以最后一个值，就会得到预期的结果。