A differential-algebraic equation is defined by an implicit functionf(du,u,p,t)=0. All of the controls are the same as theother examples, except here you define a function which returns theresiduals for each part of the equation to define the DAE. The initialvalueu0 and the initial derivativedu0 arerequired, though they do not necessarily have to satisfyf(known as inconsistent initial conditions). The methods willautomatically find consistent initial conditions. In order for this tooccur,differential_vars must be set. This vector stateswhich of the variables are differential (have a derivative term), withfalse meaning that the variable is purely algebraic.
This example shows how to solve the Robertson equation:
f<-function (du,u,p,t) { resid1=-0.04*u[1]+1e4*u[2]*u[3]- du[1] resid2=+0.04*u[1]-3e7*u[2]^2-1e4*u[2]*u[3]- du[2] resid3= u[1]+ u[2]+ u[3]-1.0c(resid1,resid2,resid3)}u0<-c(1.0,0,0)du0<-c(-0.04,0.04,0.0)tspan<-c(0.0,100000.0)differential_vars<-c(TRUE,TRUE,FALSE)prob<- de$DAEProblem(f,du0,u0,tspan,differential_vars=differential_vars)sol<- de$solve(prob)udf<-as.data.frame(t(sapply(sol$u,identity)))plotly::plot_ly(udf,x = sol$t,y =~V1,type ='scatter',mode ='lines')%>% plotly::add_trace(y =~V2)%>% plotly::add_trace(y =~V3)Additionally, an in-place JIT compiled form forf can beused to enhance the speed:
f= JuliaCall::julia_eval("function f(out,du,u,p,t) out[1] = - 0.04u[1] + 1e4*u[2]*u[3] - du[1] out[2] = + 0.04u[1] - 3e7*u[2]^2 - 1e4*u[2]*u[3] - du[2] out[3] = u[1] + u[2] + u[3] - 1.0end")u0<-c(1.0,0,0)du0<-c(-0.04,0.04,0.0)tspan<-c(0.0,100000.0)differential_vars<-c(TRUE,TRUE,FALSE)JuliaCall::julia_assign("du0", du0)JuliaCall::julia_assign("u0", u0)JuliaCall::julia_assign("p", p)JuliaCall::julia_assign("tspan", tspan)JuliaCall::julia_assign("differential_vars", differential_vars)prob= JuliaCall::julia_eval("DAEProblem(f, du0, u0, tspan, p, differential_vars=differential_vars)")sol= de$solve(prob)