如何在 Matlab 中求解常微分方程（ODE）系统

Question

我必须求解以下形式的常微分方程组：

dx/ds = 1/x * [y* (g + s/y) - a*x*f(x^2,y^2)]
dy/ds = 1/x * [-y * (b + y) * f()] - y/s - c

其中x，y是我需要找出的变量，s是自变量；其余的都是常数。到目前为止，我尝试使用 ode45 解决此问题，但没有成功：

y = ode45(@yprime, s, [1 1]);

function dyds = yprime(s,y)
    global g a v0 d
    dyds_1 = 1./y(1) .*(y(2) .* (g + s ./ y(2)) - a .* y(1) .* sqrt(y(1).^2 + (v0 + y(2)).^2));
    dyds_2 = - (y(2) .* (v0 + y(2)) .* sqrt(y(1).^2 + (v0 + y(2)).^2))./y(1) - y(2)./s - d;
   dyds = [dyds_1; dyds_2];
return

其中@yprime 具有方程组。我收到以下错误消息：

YPRIME 返回一个长度为 0 的向量，但初始的长度 条件向量为 2。 YPRIME 返回的向量和初始 条件向量必须具有相同数量的元素。

有什么想法吗？ 谢谢

Answer 1

当然，你应该看看你的函数yprime。使用一些与您的问题共享微分状态变量数量的简单模型，看看这个例子。

function dyds = yprime(s, y)
    dyds = zeros(2, 1);
    dyds(1) = y(1) + y(2);
    dyds(2) = 0.5 * y(1);
end

yprime 必须返回一个包含右侧两个值的列向量。输入参数 s 可以忽略，因为您的模型与时间无关。您展示的示例有些困难，因为它不是 dy/dt = f(t, y) 的形式。作为第一步，您必须重新排列方程式。这将有助于将x 重命名为y(1) 并将y 重命名为y(2)。

另外，您确定您的global g a v0 d 不为空吗？如果这些变量中的任何一个仍未初始化，则您将状态变量与一个空矩阵相乘，最终导致返回一个空向量 dyds。这可以用

进行测试

assert(~isempty(v0), 'v0 not initialized');

在yprime 中，或者您可以使用调试断点。

Answer 2

ODE 求解器的语法是 [s y]=ode45(@yprime, [1 10], [2 2])

并且您不需要在您的情况下进行元素操作，即使用 .* 而不是 *