MATLAB 中 ODE 解的导数值错误
Wrong value for the derivative of the solution of an ODE in MATLAB
尝试使用 here 中描述的经典 Runge-Kutta 方法求解非线性 ODE 的柯西问题(参见 ode4
)。 运行这部分之后
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。如何修正inaccuracy/mistake?非常感谢任何帮助。
函数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
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
无法计算出您似乎期望的差商。如果您阅读文档,您应该会发现它只计算差异。而第一个差异 return 只是第一个值 1.66652642851402e-05
.
导致 ode4
算法产生奇怪结果的直接错误是 eqGreen
在应该 return 行向量的地方生成列向量。由于初始值是一个行向量,在 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
缩放最后一个,您将得到预期的结果。
尝试使用 here 中描述的经典 Runge-Kutta 方法求解非线性 ODE 的柯西问题(参见 ode4
)。 运行这部分之后
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。如何修正inaccuracy/mistake?非常感谢任何帮助。
函数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
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
无法计算出您似乎期望的差商。如果您阅读文档,您应该会发现它只计算差异。而第一个差异 return 只是第一个值 1.66652642851402e-05
.
导致 ode4
算法产生奇怪结果的直接错误是 eqGreen
在应该 return 行向量的地方生成列向量。由于初始值是一个行向量,在 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
缩放最后一个,您将得到预期的结果。