将 fsolve 嵌入到 ODE 中，而不是在 Scilab 中使用 DAE 求解器

Question

我试图在 Scilab 中嵌入一个 Fsolve 来求解这个非线性系统。

我已经解决了 DAE 的问题，所以我知道会发生什么，但我正在努力嵌入 Fsolve。

这是代码的完整副本，包括 DAE。

我不确定在哪里嵌入 fsolve 函数。

//dassl approach
//given data
p=[20.086, 8100, 20.086, 20.086, 4050, 1E-17, 1E-11, 1E-17] //mol/kgh
ynames = ['y1' 'y2' 'y3' 'y4' 'y5' 'y6' 'y7' 'y8' 'y9' 'y10']
y0=[1.5776, 0.32, 0, 0, 0, 0.0131, 4.0E-06, 4.0E-06, 0, 0]'//initial conditions 
t0=0
t=linspace(0,1,1000)

dy0 = [0 0 0 0 0 0 0 0 0 0]'
function [f,r]=BatchRXorDassl(t,y,dy)
    f(1)=-p(3)*y(2)*y(8)-dy(1)
    f(2)=-p(1)*y(2)*y(6)+p(2)*y(10)-p(3)*y(2)*y(8)-dy(2)
    f(3)=p(3)*y(2)*y(8)+p(4)*y(4)*y(6)-p(5)*y(9)-dy(3)
    f(4)=-p(4)*y(4)*y(6)+p(5)*y(9)-dy(4)
    f(5)=p(1)*y(2)*y(6)-p(2)*y(10)-dy(5)
    f(6)=-p(1)*y(2)*y(6)-p(4)*y(4)*y(6)+p(2)*y(10)+p(5)*y(9)-dy(6)
    f(7)=-0.0131+y(6)+y(8)+y(9)+y(10)-y(7)
    f(8)=p(7)*y(1)-y(8)*(p(7)+y(7))
    f(9)=p(8)*y(3)-y(9)*(p(8)+y(7))
    f(10)=p(6)*y(5)-y(10)*(p(6)+y(7))
    r=0
endfunction

y=dassl([y0,dy0],t0,t,BatchRXorDassl)
y=y'
tplot = y(:,1)
yplot = y(:,2:11)

clf(11), scf(11)
plot(tplot,yplot)
xlabel('t (s)')
ylabel('C')
legend(ynames,-4)
xtitle('Batch Reactor Concentration Profiles')

////embedded fsolve approach
////given data
//p=[20.086, 8100, 20.086, 20.086, 4050, 1E-17, 1E-11, 1E-17] //mol/kgh
//ynames = ['y1' 'y2' 'y3' 'y4' 'y5' 'y6' 'y7' 'y8' 'y9' 'y10']
//y0=[1.5776, 0.32, 0, 0, 0, 0.0131, 4.0E-06, 4.0E-06, 0, 0]'//initial conditions 
//t0=0
//t=linspace(0,1,1000)
//function ff=fsolvve(y)
//    y1 = y(1)
//    y2 = y(2)
//    y3 = y(3)
//    y4 = y(4)
//    y5 = y(5)
//    y6 = y(6)
//    y7 = y(7)
//    y8 = y(8)
//    y9 = y(9)
//    y10 = y(10)
//    ff(1) = -0.0131+y(6)+y(8)+y(9)+y(10)-y(7)
//    ff(2) = p(7)*y(1)-y(8)*(p(7)+y(7))
//    ff(3) = p(8)*y(3)-y(9)*(p(8)+y(7))
//    ff(4) = p(6)*y(5)-y(10)*(p(6)+y(7))
//    ff(5) = -p(3)*y(2)*y(8)
//    ff(6) = -p(1)*y(2)*y(6)+p(2)*y(10)-p(3)*y(2)*y(8)
//    ff(7) = p(3)*y(2)*y(8)+p(4)*y(4)*y(6)-p(5)*y(9)
//    ff(8) = -p(4)*y(4)*y(6)+p(5)*y(9)
//    ff(9) = p(1)*y(2)*y(6)-p(2)*y(10)
//    ff(10) = -p(1)*y(2)*y(6)-p(4)*y(4)*y(6)+p(2)*y(10)+p(5)*y(9)
//endfunction
//yfsolve0 = [0 0 0 0 0 0 0 0 0 0]'
//yfsolve = fsolve(yfsolve0,fsolvve)

Answer 1

如果您大大放宽 ode 的默认 RTOL 和 ATOL，则以下代码可以满足您的需求：

p = [20.086, 8100, 20.086, 20.086, 4050, 1E-17, 1E-11, 1E-17] //mol/kgh
ynames = ['y1' 'y2' 'y3' 'y4' 'y5' 'y6' 'y7' 'y8' 'y9' 'y10']
y0 = [1.5776, 0.32, 0, 0, 0, 0.0131, 4.0E-06, 4.0E-06, 0, 0]'//initial condition 
t0 = 0
t = linspace(0,1,1000)

function dy = rhs(t,y1)
    y2 = fsolve([1;1;1;1],list(cstr,y1));
    y = [y1;y2];
    dy = [-p(3)*y(2)*y(8)
           p(1)*y(2)*y(6)+p(2)*y(10)-p(3)*y(2)*y(8)
           p(3)*y(2)*y(8)+p(4)*y(4)*y(6)-p(5)*y(9)
          -p(4)*y(4)*y(6)+p(5)*y(9)
           p(1)*y(2)*y(6)-p(2)*y(10)
          -p(1)*y(2)*y(6)-p(4)*y(4)*y(6)+p(2)*y(10)+p(5)*y(9)];
endfunction

function g = cstr(y2,y1)
    g = [-0.0131+y1(6)+y2(2)+y2(3)+y2(4)-y2(1)
         p(7)*y1(1)-y2(2)*(p(7)+y2(1))
         p(8)*y1(3)-y2(3)*(p(8)+y2(1))
         p(6)*y1(5)-y2(4)*(p(6)+y2(1))];
endfunction

y1 = ode(y0(1:6),t0,t,1e-3,1e-3,rhs)

必须这样做表明以这种方式进行是无稽之谈。请注意 fsolve 使用迭代算法，使用这种方法您甚至无法为其获取足够的初始向量。如果您绝对想使用 ode，我可以提出一种使用隐式微分的更好方法。

将 fsolve 嵌入到 ODE 中，而不是在 Scilab 中使用 DAE 求解器

Embedding an fsolve into an ODE as opposed to using DAE solver in Scilab

scilab