matlab - 求解一个简单的微分方程组

问题描述

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

笔记：

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

标签： matlaboctavenumerical-methodsdifferential-equations

让我们首先复制香草解决方案

% 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"]);

积分器失败，如报告的警告

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

重复积分直到关键时间前不久

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

然后导致 ODEu

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'"]);

其中列出了uas的零点0.4949319379979706, 2.886092605590324，再次确认了警告的原因，并给出了情节

matlab - 求解一个简单的微分方程组

问题描述

解决方案

推荐阅读