Movatterモバイル変換


[0]ホーム

URL:


Solving Differential-Algebraic Equations(DAE) in R with diffeqr

Chris Rackauckas

2024-12-04

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)
daes
daes

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)

[8]ページ先頭

©2009-2025 Movatter.jp