如何在Scilab中导出ODE求解时的右端函数(rhs)全时间步数据？

问题背景 #

在Scilab中编写生物行为的ODE微分求解函数时，需要捕获dy(4)和dy(5)的全时间步数据，但目前仅能获取最后一个时间步的结果。控制台disp能显示所有时间步的dy，但无法将这些数据保存为变量；尝试过全局变量、额外输出、内部调用存储函数等方法均未生效。

解决方案 #

Scilab的ODE求解器（如ode()）默认仅返回最终的时间序列和状态变量，不会自动保存微分函数内的中间计算值。以下是两种可靠的实现方式：

方法一：全局变量+微分函数内数据追加 #

这种方式直接在微分函数中把需要的数据追加到全局变量中，逻辑简单直接：

1. 主脚本初始化全局变量

// 初始化全局存储变量，用于保存每一步的额外数据
global extra_data;
extra_data = [];

// 定义初始条件、参数、时间范围（根据你的实际需求修改）
y0 = [初始值1; 初始值2; 初始值3];
cte = [参数1; 参数2; ... 参数13];
VrSTR = [你的启动子强度数据];
tspan = [0, 100]; // 模拟的时间区间

// 调用ODE求解器，传入微分函数及参数
[t, y] = ode(y0, tspan, list(diferentialSolver, cte, VrSTR));

// 查看最终存储的全时间步额外数据
disp(extra_data);

2. 修改微分函数，追加数据到全局变量

function [dy] = diferentialSolver(t,y,cte,VrSTR)
    global extra_data; // 声明全局变量，确保能访问主脚本中的存储容器
    
    // 原有计算逻辑保持不变
    [R_pfree] = cubicRoot(cte,y)
    R_pf=[]
    
    if abs(imag(R_pfree(6)))<1D-10 then
        R_pf = abs(real(R_pfree(6)))
    elseif abs(imag(R_pfree(7)))<1.D-10 then
        R_pf = abs(real(R_pfree(7)))
    elseif abs(imag(R_pfree(8)))<1.D-10
        R_pf = abs(real(R_pfree(8)))
    end

    V = y(1)
    R_T= y(2)
    if length(y)==3 then
        F = y(3)
    else
        F = y(3:$)
    end
    
    k_fR=cte(1)
    k_fB=cte(2)
    k_dR=cte(3)
    k_dB=cte(4)
    Q = cte(5)
    r_T = cte(6)
    b_T=cte(7)
    K_dQ=cte(8)
    K_dQ2=cte(9)
    K_dR=cte(10)
    K_dB=cte(11)
    r=cte(12)
    k= cte(13)
    vr=VrSTR(:)
    
    R_gfree =r_T*K_dR/(K_dR+R_pf)
    BM_gfree=b_T*K_dB/(K_dB +R_pf)
    
    dy(1)= r*V*(1-(V/k))
    dy(2)= (k_fR*vr*R_gfree)-(k_dR*R_pf)-((y(2)/V)*dy(1))
    
    if length(y)==3 then
        dy(3)= k_fB*BM_gfree - k_dB*F-F/V*dy(1)
        dy(4)=R_gfree;
        dy(5)=BM_gfree;
        // 将当前时间、dy(4)、dy(5)追加到全局变量
        extra_data = [extra_data; t, dy(4), dy(5)];
    else
        for i = 1:length(F)
            dy(2+i)= k_fB*BM_gfree - k_dB*F(i)-(F(i)/V*dy(1))
        end
        dy(10)=R_gfree;
        dy(11)=BM_gfree;
        // 多荧光变量场景下同样追加数据
        extra_data = [extra_data; t, dy(10), dy(11)];
    end
endfunction

方法二：使用ODE求解器的自定义输出函数 #

这种方式将数据存储逻辑与微分函数分离，更模块化，适合复杂场景：

1. 主脚本定义输出函数

// 自定义输出函数，每一步求解完成后被ODE求解器调用
function my_output(t, y, dy, flag)
    global extra_data;
    if flag == 0 then
        // 初始化存储容器
        extra_data = [];
    elseif flag == 2 then
        // 每一步求解完成后，提取并存储需要的数据
        // 从主脚本获取参数（避免重复传递）
        cte = get('cte', 'caller');
        VrSTR = get('VrSTR', 'caller');
        
        // 重新计算需要的额外数据（或直接从dy中提取）
        [R_pfree] = cubicRoot(cte,y);
        R_pf=[];
        if abs(imag(R_pfree(6)))<1D-10 then
            R_pf = abs(real(R_pfree(6)));
        elseif abs(imag(R_pfree(7)))<1.D-10 then
            R_pf = abs(real(R_pfree(7)));
        elseif abs(imag(R_pfree(8)))<1.D-10 then
            R_pf = abs(real(R_pfree(8)));
        end
        R_gfree = cte(6)*cte(10)/(cte(10)+R_pf);
        BM_gfree = cte(7)*cte(11)/(cte(11)+R_pf);
        
        // 追加到存储容器
        extra_data = [extra_data; t, R_gfree, BM_gfree];
    end
endfunction

// 初始化全局变量
global extra_data;
extra_data = [];

// 定义初始条件、参数等（同方法一）
y0 = [初始值1; 初始值2; 初始值3];
cte = [参数1; 参数2; ... 参数13];
VrSTR = [你的启动子强度数据];
tspan = [0, 100];

// 调用ODE求解器时指定自定义输出函数
[t, y] = ode(y0, tspan, list(diferentialSolver, cte, VrSTR), [], [], my_output);

// 查看结果
disp(extra_data);

常见问题说明 #

• 之前使用全局变量失败，大概率是因为未在主脚本和微分函数中同时声明global，或者初始化时未将变量设为空矩阵，导致数据无法正确追加。
• 微分函数的额外输出无效，是因为Scilab的ODE求解器只会接收并处理dy作为输出，其他返回值只会在最后一步被返回，无法捕获中间步骤的数据。

内容的提问来源于stack exchange，提问作者Franco Endrigo