求解一个简单的微分方程组答案

【问题标题】：Solving a simple system of differential equations求解一个简单的微分方程组
【发布时间】：2020-02-17 21:13:12
【问题描述】：

我将如何使用 Octave 数值求解以下简单的微分方程组？

注意：

我使用限定词“简单”，因为据我了解，系统是一阶且不耦合。
我已经尝试了每一个在线尝试解决此问题的方法和脚本，包括here， here 和 here。在所有选择中，我要么上吊，无响应的 Octave，提示说明“重复收敛失败”，我手动调整的建议错误初始步长和最大步长（我确实尝试过，但没有有用），或者最初看起来像解决方案的东西，因为没有错误，但绘制解决方案会显示一个空白图表
Octave 提供了等效的 Matlab 例程，我尝试了各种例程ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i，当然还有, Octaves 自己的lsode 命令，都给出了上述相同的错误。

【问题讨论】：

标签： matlab octave numerical-methods differential-equations

【解决方案1】：

让我们首先复制原版解决方案

% z = [x,y]
f = @(t,z) [ z(1).^2+t; z(1).*z(2)-2 ];
z0 = [ 2; 1];

[ T, Z ] = ode45(f, [0, 10], z0);

plot(T,Z); legend(["x";"y"]);

集成器失败并显示警告

警告：求解不成功。在到达tend = 10.000000 的端点之前，迭代集成循环在时间t = 0.494898 退出。如果步长变得太小，可能会发生这种情况。尝试使用命令'odeset' 减小'InitialStep' 和/或'MaxStep' 的值。

在关键时间前不久重复集成

opt = odeset('MaxStep',0.01);
[ T, Z ] = ode45(f, [0, 0.49], z0, opt);
clf; plot(T,Z); legend(["x";"y"]);

图表中的结果

我们可以看到，第一个方程中的二次项会导致增长失控。由于某种原因，求解器只识别不断减小的步长，而不识别解的失控值。

确实，第一个是 Riccati 方程，已知它在有限时间具有极点。使用典型的参数化x(t)=-u'(t)/u(t) 具有乘积/商规则的导数

x' = -u''(t)/u(t) - u'(t)* (-u'(t)/u(t)^2) = -u''(t)/u(t) + x(t)^2

然后导致 u 的 ODE

u''(t)+t*u(t)=0, u(0)=-1, u'(0)=x(0)=2,

这是一个艾里方程，带有t>0 的振荡分支。 u 的第一个根是x 的极点，没有办法将解扩展到这一点之外。

g=@(t,u) [u(2); -t.*u(1)]
u0 = [ 1; -2];

function [val,term, dir] = event(t,u)
   val = u(1);
   term = 0;
   dir = 0;
end
opt = odeset('MaxStep',0.1, 'Events', @(t,u) event(t,u));
[T,U,Te,Ue,Ie] = ode45(g,[0,4],u0,opt);
disp(Te)
clf; plot(T,U); legend(["u";"u'"]);

将u 的零列为0.4949319379979706, 2.886092605590324，再次确认警告的原因，并给出情节

【讨论】：

当你开始参数化解决方案时我失去了你（由于我自己缺乏数学理解），但如果我理解你答案的本质，本质上存在一个我们必须解决的不连续性.在那种情况下，我是否理解您是否正确，点 0.4949 之后有解决方案。例如，在 x(1) 和 y(1) 处有一个值？
不，解决方案在奇点处结束，没有有意义的延续。当然 $u$ 无处不在，所以分数 $-u'/u$ 存在于极点的两侧，但是这两个线段不相连，从奇点的左边到奇点的右边没有无穷小的步长。