与Scilab中使用DAE求解器相反，将fsolve嵌入到ODE中

如何解决与Scilab中使用DAE求解器相反，将fsolve嵌入到ODE中

我正在尝试在Scilab中嵌入Fsolve来解决此非线性系统。

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

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

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

//dassl approach
//given data
p=[20.086,8100,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.0131,4.0E-06,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,1E-17] //mol/kgh
//ynames = ['y1' 'y2' 'y3' 'y4' 'y5' 'y6' 'y7' 'y8' 'y9' 'y10']
//y0=[1.5776,0]'//initial conditions 
//t0=0
//t=linspace(0,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)

解决方法

如果您大大放松了ode的默认RTOL和ATOL，以下代码将满足您的要求：

p = [20.086,8100,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.0131,4.0E-06,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,rhs)

必须这样做表明这样做是没有意义的。请意识到fsolve使用了迭代算法，并且使用这种方法，您甚至都无法为其获取足够的初始向量。如果您绝对想使用ode，我可以提出一种更好的使用隐式微分的方法。