IDASolve 将 ODE 从 Matlab 转换为 Octave 时出错

Question

我目前正在尝试将一些脚本从 Matlab 转换为 Octave，包括 ODE 求解器。在 Matlab 上，脚本运行顺利，但 Octave 返回：

[IDA ERROR]  IDASolve
  At t = 0.00101 and h = 3.8147e-12, the corrector convergence failed repeatedly or with |h| = hmin.

为了简单起见，我只会发布相关的函数，变量在上游进一步定义：

% Create a function handle with constant parameters
DissIonWithConstants = @(xx, ab) DissIon(xx, ab, [B_d_f, k_b, j, s_d_f, E_D, R, B_i_f, s_i_f, E_I, Ma_dgl, M_mol_H, rho_d, m_H, t_el, rho_dgl, alpha0, delta0]);

# ode-solver settings
OPTs = odeset('reltol', 1e-12, 'abstol', 1e-16, 'normcontrol', 'on', 'refine', 5);

% Solve the ODE using the function handle
[xx, ab] = ode15s(DissIonWithConstants, x_span, [alpha0, delta0],OPTs);


function dab = DissIon(xx, ab, werte)
    % dab is output
    dab = zeros(2, 1);

    B_d_f = werte(1);
    k = werte(2);
    T = werte(3);
    s_d_f = werte(4);
    E_D = werte(5);
    R_universal = werte(6);
    B_i_f = werte(7);
    s_i_f = werte(8);
    E_I = werte(9);
    Ma = werte(10);
    M_mol_H = werte(11);
    rho_d = werte(12);
    m_H = werte(13);
    T_e = werte(14);
    rho = werte(15);
    alpha_start = werte(16);
    delta_start = werte(17);

    % Ionisation 
    dab(1) = ab(1) * B_i_f * T_e^(s_i_f) * exp(-E_I/(k*T_e)) * ...
             (rho*(ab(2)-ab(1))/M_mol_H - ...
             (ab(2)-ab(1))/(m_H*(1-ab(1))) * ab(1)^2 * rho^2 / M_mol_H * 2.4147 * 10^(-21) * T_e^(-3/2) * exp(E_I/(k*T_e))) * ...
             1 / (Ma * ((R_universal/M_mol_H*T*(1+alpha_start*T_e/T)/(2-delta_start))^(1/2)));

    % Dissociation
    dab(2) = (ab(2)-2*ab(1)+1) * B_d_f * T^(s_d_f) * exp(-E_D/(k*T)) * ...
             (rho*(1-ab(2))/(2*M_mol_H) - ab(2)^2 * rho^2 / (2*rho_d*M_mol_H) * exp(E_D/(k*T))) * ...
             1 / (Ma * ((R_universal/M_mol_H*T*(1+alpha_start*T_e/T)/(2-delta_start))^(1/2)));
end

我使用了 odeset，但无法运行它。

知道为什么这在 Octave 中不起作用吗？我尝试使用 GPT 进行优化，但没有成功。

任何帮助将不胜感激，谢谢！

Answer 1

您的代码在文档/自文档代码方面还有很多不足之处，但即便如此，这个代码片段还是有几件事有点“奇怪”，即使人们不需要知道它想要做什么做。我怀疑发生了以下两种情况之一：要么您的工作空间中有尚未清除的东西，要么您在没有意识到的情况下使用旧东西。或者，您以错误的方式/错误的参数调用该函数。

一方面，即使我编写变量来让代码运行，它也不会，因为你的函数在它甚至还没有定义的时候被调用。请注意，matlab 期望“脚本”函数最后出现，但 Octave 期望它们在使用之前定义。因此，每次运行此脚本时，您都在运行“旧”DissIon 函数。

在缺乏至少有关您传递的变量的“类型”和维度的任何信息的情况下，这里没什么可说的了。 但是，您可能会发现以下技术有助于帮助您调试问题。

1. 使用

   debug_on_error

   /

   debug_on_warning

   函数帮助您在生成错误/警告时停止并检查当前运行函数的状态。

2. 像这样高度数学/方程驱动的代码非常难以调试。尝试将其拆分为具有有意义名称的逻辑组件，和/或在涉及关系时使用函数（换句话说，尝试使用函数式编程方法而不是状态驱动方法，因为状态驱动方法可能会引发错误）。

3. 不要像处理“werte”那样传递需要解构的东西，而是直接传递一个结构体。这样，如果您不小心弄乱了关键字的位置，也不会导致逻辑错误。

即便如此，如果你“确实”需要八度解构赋值，八度通过逗号分隔列表支持这一点，你不必像这样手动一一完成。例如。只需收集元胞数组中的所有右侧标记（例如

C

），并按如下方式分配它：[a,b,c] = C{:}。 （有关详细信息，请参阅 https://docs.octave.org/latest/Comma_002dSeparated-Lists.html，以及

deal

命令的文档）；

最后，请记住，虽然某些函数具有相同的名称，但 Octave 中的实现可能与 Matlab 不同，存在细微的差异。因此，值得检查每种情况的文档，以确认它们是否采用相同的参数、相同的顺序、形式，并对输入做出相同的假设。