0

我在 MATLAB 中使用 ode45 运行一组 ODE，我需要保存其中一个变量（这不是导数）以供以后使用。我正在使用函数“assignin”在基础工作区中分配一个临时变量，并在每一步更新它。这似乎可行，但是，数组的大小与从 ode45 获取的解向量的大小不匹配。例如，我有以下嵌套函数：

function [Z,Y] = droplet_momentum(theta,K,G,P,zspan,Y0)

options = odeset('RelTol',1e-7,'AbsTol',1e-7);
[Z,Y] = ode45(@momentum,zspan,Y0,options);

function DY = momentum(z,y)

    DY = zeros(4,1);

    %Entrained Total Velocity
    Ve = sin(theta)*(y(4));

    %Total Relative Velocity
    Urs = sqrt((y(1) - y(4))^2 + (y(2) - Ve*cos(theta))^2 + (y(3))^2);

    %Coefficients
    PSI = K*Urs/y(1);
    PHI = P*Urs/y(1);

    %Liquid Axial Velocity
    DY(1) = PSI*sign(y(1) - y(4))*(1 + (1/6)*(abs(y(1) - y(4))*G)^(2/3));

    %Liquid Radial Velocity
    DY(2) = PSI*sign(y(2) - Ve*cos(theta))*(1 + (1/6)*(abs(y(2) - ...
        Ve*cos(theta))*G)^(2/3));

    %Liquid Tangential Velocity
    DY(3) = PSI*sign(y(3))*(1 + (1/6)*(abs(y(3))*G)^(2/3));

    %Gaseous Axial Velocity
    DY(4) = (1/z/y(4))*((PHI/z)*sign(y(1) - y(4))*(1 + ...
        (1/6)*(abs(y(1) - y(4))*G)^(2/3)) + Ve*Ve - y(4)*y(4));

    assignin('base','Ve_step',Ve);
    evalin('base','Ve_out(end+1) = Ve_step');
end

end

在上面的代码中，theta（弧度）、K（负值）、P 和 G 是常数，为了这个例子，可以取任何值。Zspan 只是 ODE 求解器的积分时间步长，Y0 是初始条件向量 (4x1)。同样，为了这个例子，这些可以取任何合理的值。现在在主文件中，使用以下内容调用该函数：

Ve_out = 0;
[Z,Y] = droplet_momentum(theta,K,G,P,zspan,Y0);
Ve_out = Ve_out(2:end);

这个方法没有被 MATLAB 抱怨，但问题是 Ve_out 的大小与 Z 或 Y 的大小不一样。原因是因为 MATLAB 对其算法多次调用 ODE 函数，所以解决方案是将比 Ve_out 略小。正如 am304 所建议的那样，我可以通过给 ode 函数一个 Z 和 Y 向量（例如 DY = 动量（Z，Y））来简单地计算 DY，但是，我需要使用“assignin”（或类似方法）来解决这个问题，因为另一个这个问题的版本在 DY 和 Ve 之间存在隐式依赖性，并且在每次迭代时计算 DY 的计算成本太高（我将在多次迭代中运行这个问题）。

4

2 回答 2

3

好的，让我们从一个 SSCCE 的快速示例开始：

function [Z,Y] = khan

options = odeset('RelTol',1e-7,'AbsTol',1e-7);
[Z,Y] = ode45(@momentum,[0 12],[0 0],options);
end

function Dy = momentum(z,y)
Dy = [0 0]';

Dy(1) = 3*y(1) + 2* y(2) - 2;
Dy(2) = y(1) - y(2);

Ve = Dy(1)+ y(2);

    assignin('base','Ve_step',Ve);
    evalin('base','Ve_out(end+1) = Ve_step;');

    assignin('base','T_step',z);
    evalin('base','T_out(end+1) = T_step;');
end

通过[Z,Y] = khan作为命令行运行，我得到了一个完整的功能代码来演示您的问题，而没有相关的所有令人头疼的问题。我对此的耐心已经耗尽：生活和学习。

这似乎可行，但是，数组的大小与从 ode45 获取的解向量的大小不匹配

请注意，我在提取时间变量的代码中添加了两行。在命令提示符下，只需运行以下命令即可了解发生了什么：

Ve_out = [];
T_out = [];
[Z,Y] = khan;
size (Z)
size (T_out)
size (Ve_out)
plot (diff(T_out))

ans =
   109     1

ans =   
     1   163

ans =
     1   163

在此处输入图像描述

基本上 ode45 是一种迭代算法，这意味着它会定期正确（这就是您经常看到 diff(T) = 0 的原因）。你不能强迫算法做你想做的事，你必须忍受它。

所以你的选择是
1.使用固定步长算法
2. 在 ode45 算法完成其肮脏的工作后，有一个函数调用可以重现你想要的东西。（am304 的解决方案）
3. 用时间函数收集数据，然后有一个算法解析所有内容以去除多余的数据。

于 2013-05-22T18:47:47.913 回答
1

你不能这样做吗？显然检查矩阵/向量的大小是否正确并相应地修改代码。

[Z,Y] = droplet_momentum2(theta,K,G,P,zspan,Y0);
DY = momentum(Z,Y);
Ve = sin(theta)*(0.5*z*DY(4) + y(4));

即，一旦求解了 ODE，就计算导数DY作为 和 的函数Z（Y刚刚由 ODE 求解），最后计算Ve。

于 2013-05-22T07:49:16.153 回答