除了解值之外，还有没有办法导出解中的 ODE rhs 值

Question

我正在编写描述生物行为的代码，除了 ODE 已经给出的输出之外，我还需要数据，这在下面更好地解释：

function [dy] = diferentialSolver(t,y,cte,VrSTR)
    //[R_pfree]= newtonBissection(cte,y)
    
    [R_pfree] =(cubicRoot(cte,y))
    [cData] = [R_pfree]
    //disp(cubicRoot(cte,y))
    R_pf=[]
    
    // this will select the root in that is a real number, meaning that will always have two complex roots in this system
    if abs(imag(R_pfree(6)))<1D-10 then
        R_pf = abs(real(R_pfree(6)))
        disp("here1")
    elseif (abs(imag(R_pfree(7))))<1.D-10 then
        R_pf = abs(real(R_pfree(7)))
        disp("here2")
    elseif abs(imag(R_pfree(8)))<1.D-10
        R_pf = abs(real(R_pfree(8)))
        disp("here3")
    end

    //disp( newtonBissection(cte,y))
    //redundant to simplify programming
    //[diffdata] = return [R_pfree]
    //R_pf= R_pfree
    V = y(1)
    R_T= y(2)
    // this if statement is to make sure the y(3) can be call either as a matrix or as a single variable.
    if length(y)==3 then
        F = y(3)
    else
        [F]=y(3:$)
    end
    //cte=abs(cte)
    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) //mi(v)
    k= cte(13)
    vr=VrSTR(:) //promoter strength 
    

    //IGR= 2 Intrinsic Growth Rate, still testing this one.
    //Ext= 1.2    //extinction rate. Not on use yet.
    
    R_gfree =r_T*K_dR/(K_dR+R_pf)
    //disp(R_gfree)    //Free repressor gene concentration at time t
    BM_gfree=b_T*K_dB/(K_dB +R_pf)
    
    //disp(BM_gfree)
   // disp(list(["Rpfree: " + string(R_pfree)],["Rgfree: " + string(R_gfree)], ["BM_gfree: " + string(BM_gfree)]))
    //dy(1)= V*%e^r*(1-V/k)        // Volume over time eq Ricker Model
    
    //the differential equations begins here
    dy(1)= r*V*(1-(V/k))    // Volume over time eq
    
    //does simulate the stabilization of the system well
     dy(2)= (k_fR*vr*R_gfree)-(k_dR*R_pf)-((y(2)/V)*dy(1))    //total repressor over time eq
    // this if statement is to make sure the y(3) can be call either as a matrix or
    //as a single variable.
    if length(y)==3 then
        dy(3)= k_fB*BM_gfree - k_dB*F-F/V*dy(1)    //fluorescence over time
        dy(4)=r_T*K_dR/(K_dR+R_pf)
        dy(5)=b_T*K_dB/(K_dB +R_pf)
    else
        for i = 1:length(F)
            //here we call F as an vector, since we accept that F is a submatrix of y(),
            //containing multiple values of Fluorescence.
            dy(2+i)= k_fB*BM_gfree - k_dB*F(i)-(F(i)/V*dy(1))    //fluorescence over time
            dy(10)=r_T*K_dR/(K_dR+R_pf)
            dy(11)=b_T*K_dB/(K_dB +R_pf)
        end
    end

    disp(dy);
endfunction

我正在尝试接收 dy 的 4 和 5，但我只输出最后一个时间步，我想要所有步骤，这可能吗？ （我的意思是应该是，因为 disp 实际上在控制台中及时显示了所有 dy，但 scilab 显然无法保存为变量）。

如果有人可以帮助我，我将非常感激！

我尝试了全局变量，只及时保存了最后一步。 我也试过上面这个方法，甚至没有保存其他dy的。 我尝试在函数中添加另一个输出，如“function [dy, ProteinFree] = differentialSolver(t,y,cte,VrSTR)”，但没有按预期工作。 另外，我尝试通过在 differentialSolver 函数内部调用它来保存 ODE 不应该干扰的另一个函数，但它不起作用。 现在我已经没有想法了。

Answer 1

在 Scilab (2023.1) 的实际版本中，产生 ode + 解时的实际 rhs 值的最快方法是将其重新表示为带有

dassl

的隐式方程（请参阅 https://help.scilab）。 org/dassl）。这是一个简单线性颂的示例：

// original ODE rhs
function yd = rhs(t,y)
  yd = [0 1;-1 0]*y;
end

// equivalent DAE residual
function [r,flag] = res(t,y,yd)
  flag = 0;
  r = yd-rhs(t,y);
end

t = linspace(0,5,100);

// solve with ode function
y0 = [0;1];
y = ode(y0,0,t,rhs);

// solve with dassl and extract y, dy
sol = dassl(y0,0,t,res);
y = sol(2:3,:);
dy = sol(4:5,:);

您还可以在

rhs

调用之后调用

ode

，如果

rhs

无法对多个

y

向量的计算进行向量化（这是您的特定 rhs 的情况），这将产生额外的成本：

y = ode(y0,0,t,rhs);
dy = zeros(y);
for i = 1:size(y,2)
  dy(:,i) = rhs(t(i),y(:,i));
end