Tutorials for doing scientific machine learning (SciML) and high-performance differential equation solving with open source software.

  • By SciML Open Source Scientific Machine Learning
  • Last update: Dec 16, 2022
  • Comments: 16

SciMLTutorials.jl: Tutorials for Scientific Machine Learning and Differential Equations

Build status

Join the chat at https://gitter.im/JuliaDiffEq/Lobby

SciMLTutorials.jl holds PDFs, webpages, and interactive Jupyter notebooks showing how to utilize the software in the SciML Scientific Machine Learning ecosystem. This set of tutorials was made to complement the documentation and the devdocs by providing practical examples of the concepts. For more details, please consult the docs.

Interactive Notebooks

To run the tutorials interactively via Jupyter notebooks, install the package and open the tutorials like:

using Pkg
pkg"add https://github.com/SciML/SciMLTutorials.jl"
using SciMLTutorials
SciMLTutorials.open_notebooks()

Video Tutorial

Video Tutorial

Table of Contents

Contributing

First of all, make sure that your current directory is SciMLTutorials. All of the files are generated from the Weave.jl files in the tutorials folder. To run the generation process, do for example:

using Pkg, SciMLTutorials
cd(joinpath(dirname(pathof(SciMLTutorials)), ".."))
Pkg.pkg"activate ."
Pkg.pkg"instantiate"
SciMLTutorials.weave_file("introduction","01-ode_introduction.jmd")

To generate all of the notebooks, do:

SciMLTutorials.weave_all()

If you add new tutorials which require new packages, simply updating your local environment will change the project and manifest files. When this occurs, the updated environment files should be included in the PR.

Github

https://github.com/SciML/SciMLTutorials.jl

Comments(16)

  • 1

    problem with animation

    Sorry for bothering you with this but I am trying to put together a notebook that animates the solution to the diffusion equation. The code below works well for plotting slices but the animation seems to go on and on and never make an animation.

    Any suggestions as to what I'm doing wrong?

    Needless to say, I am happy to share the notebook in case anyone might find it helpful.

    # Import packages
    using DiffEqOperators, DifferentialEquations, Plots
    
    x = collect(-pi : 2pi/511 : pi);
    u0 = -(x - 0.5).^2 + 1/12;
    A = DerivativeOperator{Float64}(2,2,2pi/511,512,:Dirichlet,:Dirichlet;BC=(u0[1],u0[end]));
    
    prob1 = ODEProblem(A, u0, (0.,10.));
    sol1 = solve(prob1, dense=false, tstops=0:0.01:0.1);
        
    # try to plot the solution at different time points using
    plot(x, [sol1(i) for i in 0:1:10])
    
    animate(sol1,"tmp.gif")
    
    
  • 2

    Enumerate tutorials and rebuild

    This PR is the continuation of #61. In a first step, I have enumerated all the tutorial files in tutorials/, html/, notebook/, pdf/, script/ as suggested in #60. The second commit rebuilds all the tutorials that could be built without unintended Weave.jl warnings/errors.

    The following tutorials could not be built without errors, so I will add them in seperate commits to this PR, and we can figure out what to do:

    • introduction/05-formatting_plots
    • models/03-diffeqbio_I_introduction
    • models/04-diffeqbio_II_networkproperties
    • models/06-pendulum_bayesian_interface
    • models/07-outer_solar_system
    • ode_extras/ModelingToolkit
    • type_handling/01-number_types
    • type_handling/03-unitful
  • 3

    Plot solutions in notebook about number with uncertainties

    The latest version of Measurements.jl supports Plots.jl and I took the occasion to update the related tutorial to show this feature (actually, I decided to implement this feature exactly for this tutorial).

    Only problem is that it doesn't seem to be possible to simply do plot(sol), I get some errors. Any idea why?

  • 4

    Example Radioactive Decay (Numbers with uncertainties) not working in Julia 1.0

    Source: http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqTutorials.jl/blob/master/PhysicalModels/NumberUncertainties.ipynb

    If I execute

    # Half-life of radiocarbon, in thousands of years
    c = 5.730 ± 0.040
    
    #Setup
    u₀ = 1 ± 0
    tspan = (0 ± 0, 1 ± 0)
    
    #Define the problem
    radioactivedecay(u,p,t) = -c * u
    
    #Pass to solver
    prob = ODEProblem(radioactivedecay, u₀, tspan)
    sol = solve(prob, Tsit5())
    

    in Julia 1.0 I get:

    MethodError: Measurement{Float64}(::Rational{Int64}) is ambiguous. Candidates:
      (::Type{Measurement{T}})(x::S) where {T, S} in Measurements at C:\Users\carsten\.julia\packages\Measurements\57KtG\src\Measurements.jl:61
      (::Type{T})(x::Rational{S}) where {S, T<:AbstractFloat} in Base at rational.jl:90
    Possible fix, define
      (::Type{Measurement{T}})(::Rational{S})
    
    Stacktrace:
     [1] #__init#203(::Int64, ::Array{Measurement{Float64},1}, ::Array{Measurement{Float64},1}, ::Array{Measurement{Float64},1}, ::Nothing, ::Bool, ::Nothing, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Measurement{Float64}, ::Bool, ::Rational{Int64}, ::Nothing, ::Nothing, ::Int64, ::Rational{Int64}, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Nothing, ::Nothing, ::Int64, ::Measurement{Float64}, ::Float64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::typeof(LinearAlgebra.opnorm), ::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), ::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), ::Nothing, ::Bool, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Measurement{Float64},Tuple{Measurement{Float64},Measurement{Float64}},false,Nothing,ODEFunction{false,typeof(radioactivedecay),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at C:\Users\carsten\.julia\packages\OrdinaryDiffEq\nZPqS\src\solve.jl:80
     [2] __init(::ODEProblem{Measurement{Float64},Tuple{Measurement{Float64},Measurement{Float64}},false,Nothing,ODEFunction{false,typeof(radioactivedecay),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at C:\Users\carsten\.julia\packages\OrdinaryDiffEq\nZPqS\src\solve.jl:61
     [3] #__solve#202(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::ODEProblem{Measurement{Float64},Tuple{Measurement{Float64},Measurement{Float64}},false,Nothing,ODEFunction{false,typeof(radioactivedecay),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at C:\Users\carsten\.julia\packages\OrdinaryDiffEq\nZPqS\src\solve.jl:6
     [4] __solve(::ODEProblem{Measurement{Float64},Tuple{Measurement{Float64},Measurement{Float64}},false,Nothing,ODEFunction{false,typeof(radioactivedecay),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at C:\Users\carsten\.julia\packages\OrdinaryDiffEq\nZPqS\src\solve.jl:6 (repeats 5 times)
     [5] #solve#429(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::ODEProblem{Measurement{Float64},Tuple{Measurement{Float64},Measurement{Float64}},false,Nothing,ODEFunction{false,typeof(radioactivedecay),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5) at C:\Users\carsten\.julia\packages\DiffEqBase\b3iNk\src\solve.jl:29
     [6] solve(::ODEProblem{Measurement{Float64},Tuple{Measurement{Float64},Measurement{Float64}},false,Nothing,ODEFunction{false,typeof(radioactivedecay),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5) at C:\Users\carsten\.julia\packages\DiffEqBase\b3iNk\src\solve.jl:21
     [7] top-level scope at In[6]:13
    
  • 5

    Tutorial - Heat Equation Finite Element Method

    Hi,

    the code in the heat equation tutorial gives the following error for me:

    using DifferentialEquations
    
    f(t,x,u)  = ones(size(x,1)) - .5u
    u0_func(x) = zeros(size(x,1))
    
    tspan = (0.0,1.0)
    dx = 1//2^(3)
    dt = 1//2^(7)
    
    mesh = parabolic_squaremesh([0 1 0 1],dx,dt,tspan,:neumann)
    u0 = u0_func(mesh.node)
    prob = HeatProblem(u0,f,mesh)
    
    sol = solve(prob,FEMDiffEqHeatImplicitEuler())
    
    Error displaying DiffEqPDEBase.FEMSolution{Float64,2,Array{Array{Float64,1},1},Void,Dict{Symbol,Float64},Array{Float64,1},DiffEqPDEBase.HeatProblem{false,false,DiffEqPDEBase.FEMMesh{Array{Rational{Int64},2},Array{Float64,1},Rational{Int64},Tuple{Float64,Float64}},#f,Void,DiffEqPDEBase.##21#30,DiffEqPDEBase.##22#31,Array{Float64,1},DiffEqPDEBase.##14#23,Void,Float64}}: type FEMSolution has no field retcode
    show at solution_interface.jl:34 [inlined]
    (::LastMain.Juno.##1#2{DiffEqPDEBase.FEMSolution{Float64,2,Array{Array{Float64,1},1},Void,Dict{Symbol,Float64},Array{Float64,1},DiffEqPDEBase.HeatProblem{false,false,DiffEqPDEBase.FEMMesh{Array{Rational{Int64},2},Array{Float64,1},Rational{Int64},Tuple{Float64,Float64}},#f,Void,DiffEqPDEBase.##21#30,DiffEqPDEBase.##22#31,Array{Float64,1},DiffEqPDEBase.##14#23,Void,Float64}}})(::Base.AbstractIOBuffer{Array{UInt8,1}}) at types.jl:11
    #sprint#228(::Void, ::Function, ::Int64, ::Function) at io.jl:66
    Type at types.jl:51 [inlined]
    Type at types.jl:52 [inlined]
    render(::LastMain.Juno.Editor, ::DiffEqPDEBase.FEMSolution{Float64,2,Array{Array{Float64,1},1},Void,Dict{Symbol,Float64},Array{Float64,1},DiffEqPDEBase.HeatProblem{false,false,DiffEqPDEBase.FEMMesh{Array{Rational{Int64},2},Array{Float64,1},Rational{Int64},Tuple{Float64,Float64}},#f,Void,DiffEqPDEBase.##21#30,DiffEqPDEBase.##22#31,Array{Float64,1},DiffEqPDEBase.##14#23,Void,Float64}}) at types.jl:13
    render′(::LastMain.Juno.Editor, ::DiffEqPDEBase.FEMSolution{Float64,2,Array{Array{Float64,1},1},Void,Dict{Symbol,Float64},Array{Float64,1},DiffEqPDEBase.HeatProblem{false,false,DiffEqPDEBase.FEMMesh{Array{Rational{Int64},2},Array{Float64,1},Rational{Int64},Tuple{Float64,Float64}},#f,Void,DiffEqPDEBase.##21#30,DiffEqPDEBase.##22#31,Array{Float64,1},DiffEqPDEBase.##14#23,Void,Float64}}) at errors.jl:106
    (::LastMain.Atom.##50#54{String})() at eval.jl:57
    macro expansion at eval.jl:53 [inlined]
    (::LastMain.Atom.##48#52{Dict{String,Any}})() at task.jl:80
    

    Running Julia v0.6.1 with DifferentialEquations v3.0.2

  • 6

    Add more dependencies to DiffEqTutorials?

    On a fresh julia-1.1.0 install, I followed the "Getting started guide" for DifferentialEquations:

    using Pkg
    Pkg.add("DifferentialEquations")
    using DifferentialEquations
    using Pkg
    pkg"add https://github.com/JuliaDiffEq/DiffEqTutorials.jl"
    using DiffEqTutorials
    DiffEqTutorials.open_notebooks()
    

    So far I've needed to add IJulia, MatrixFactorizations, DiffEqParamEstim, Plots, Optim ... and perhaps a few more to come (haven't gone through the tutorial yet).

    I am somewhat perplexed since I see these requirements listed here: https://github.com/JuliaDiffEq/DiffEqTutorials.jl/blob/master/Project.toml

    Is that pkg"add ..." special in the sense it doesn't pull the dependencies?

  • 7

    Update most if not all notebooks

    @ChrisRackauckas Can you or someone please update most if not all notebooks in DiffEqTutorials to run smoothly on current Julia v1.0.1 and DifferentialEquations v5.3.1. For a fruitful dissemination of Julia DiffEq, it is important that this collection of notebooks runs flawlessly.

  • 8

    DiffEqBiological packages outdated

    The tutorials https://tutorials.sciml.ai/html/models/03-diffeqbio_I_introduction.html and following are outdated. They still refer to DiffEqBiological instead of Catalyst. On my system, I started to fix the tutorials and I came across several issues, some of which I could fix. I will try to create a comprehensive list of errors if you are interested in keeping these tutorials. In my opinion this would make sense, because the examples in the Catalyst docs are by far not as detailed.

  • 9

    updated kepler problem tutorial

    1. Corrected arguments order in calls of qdot and pdot.
    2. There are different orders of variables u and v in ODEProblem and DynamicalODEProblem. Moreover, I guess, the order was changed at some time of these problems' existence. I have updated the plotting functions with the current state of things.
    3. Small remark: it seems like in the documentation a wrong order of u and v variables noted in the declaration of f1 and f2.
  • 10

    ExtraFEMFeatures/Heat Equation

    I am very new to Julia and going through some of the neat tutorials people have put together.

    I installed Julia 0.6.0 on Ubunbu by downloading the generic version, as requested.

    When I run the Heat Equation example the first line seems to go okay but the second has a MethodError, suggesting a method might be too new.

    To be able to use this example should I use the most recent (bleeding edge) version? I am happy to do that but wanted to check to see if this is what is necessary.

    MethodError: no method matching setcharheight(::Float64) The applicable method may be too new: running in world age 24141, while current world is 24142. Closest candidates are: setcharheight(::Real) at /home/fpoulin/.julia/v0.6/GR/src/GR.jl:1010 (method too new to be called from this world context.)

    Stacktrace:
     [1] #gr_set_font#316(::Symbol, ::Symbol, ::ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}}, ::Int64, ::Function, ::Plots.Font) at /home/fpoulin/.julia/v0.6/Plots/src/backends/gr.jl:374
     [2] (::Plots.#kw##gr_set_font)(::Array{Any,1}, ::Plots.#gr_set_font, ::Plots.Font) at ./<missing>:0
     [3] gr_set_xticks_font(::Plots.Subplot{Plots.GRBackend}) at /home/fpoulin/.julia/v0.6/Plots/src/backends/gr.jl:561
     [4] _update_min_padding!(::Plots.Subplot{Plots.GRBackend}) at /home/fpoulin/.julia/v0.6/Plots/src/backends/gr.jl:612
     [5] _collect(::Array{RecipesBase.AbstractLayout,2}, ::Base.Generator{Array{RecipesBase.AbstractLayout,2},Plots.#_update_min_padding!}, ::Base.EltypeUnknown, ::Base.HasShape) at ./array.jl:454
     [6] _update_min_padding!(::Plots.GridLayout) at /home/fpoulin/.julia/v0.6/Plots/src/layouts.jl:310
     [7] prepare_output(::Plots.Plot{Plots.GRBackend}) at /home/fpoulin/.julia/v0.6/Plots/src/plot.jl:258
     [8] show(::IOStream, ::MIME{Symbol("image/png")}, ::Plots.Plot{Plots.GRBackend}) at /home/fpoulin/.julia/v0.6/Plots/src/output.jl:207
     [9] png(::Plots.Plot{Plots.GRBackend}, ::String) at /home/fpoulin/.julia/v0.6/Plots/src/output.jl:8
     [10] frame(::Plots.Animation, ::Plots.Plot{Plots.GRBackend}) at /home/fpoulin/.julia/v0.6/Plots/src/animation.jl:20
     [11] #animate#252(::Array{Any,1}, ::Function, ::Plots.FrameIterator, ::String) at /home/fpoulin/.julia/v0.6/Plots/src/animation.jl:43
     [12] (::RecipesBase.#kw##animate)(::Array{Any,1}, ::RecipesBase.#animate, ::Plots.FrameIterator, ::String) at ./<missing>:0 (repeats 2 times)
     [13] include_string(::String, ::String) at ./loading.jl:515
    
  • 11

    Rebuild tutorials automatically

    First draft, I expect some things to be broken.

    Here's what's supposed to happen:

    • Every 3 days, GitHub Actions fires and triggers GitLab CI
    • When GitLab CI is triggered, the setup job chooses a random folder and file, then emits a job config for that tutorial.
    • The rebuild GitLab CI job reads this generated job config and executes it, by:
      • Building the tutorial
      • If there are changes, commit and push them to a new branch
    • GitHub Actions receives notification of the new push, and opens a PR for that branch
    • You can also manually trigger a rebuild of a specific tutorial by making an issue comment with the contents !rebuild <folder>/<file>, in which case the same process happens but without the random selection
  • 12

    SciML Tutorial attempts to simultaneously use keywords vars and idxs in Plot, which is not supported.

    When I try to plot the first 250 trajectories, according to the SciML tutorial:

    https://github.com/SciML/SciMLTutorials.jl/blob/e8c8c3e51703fd08b874abbbc7b52ed895d0cb6f/tutorials/DiffEqUncertainty/01-expectation_introduction.jmd#L48

    I get this error:

    Simultaneously using keywords vars and idxs is not supported. Please only use idxs.
    
    Stacktrace:
     [1] error(s::String)
       @ Base ./error.jl:35
     [2] macro expansion
       @ ~/.julia/packages/SciMLBase/UY87O/src/solutions/solution_interface.jl:189 [inlined]
     [3] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, sol::SciMLBase.AbstractTimeseriesSolution)
       @ SciMLBase ~/.julia/packages/RecipesBase/qpxEX/src/RecipesBase.jl:289
     [4] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
       @ RecipesPipeline ~/.julia/packages/RecipesPipeline/OXGmH/src/user_recipe.jl:36
     [5] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
       @ RecipesPipeline ~/.julia/packages/RecipesPipeline/OXGmH/src/RecipesPipeline.jl:70
     [6] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
       @ Plots ~/.julia/packages/Plots/lW9ll/src/plot.jl:209
     [7] #plot#145
       @ ~/.julia/packages/Plots/lW9ll/src/plot.jl:91 [inlined]
     [8] top-level scope
       @ ~/Projects/Proposals/ASKEM/SciMLTutorials.jl/tutorials/DiffEqUncertainty/01-expectation_introduction.jmd.ipynb:2
    
  • 13

    Double Pendulum example lacks explanations

    The double pendulum example has a certain mismatch between the code and the explanations surrounding it: https://tutorials.sciml.ai/html/models/01-classical_physics.html

    It firstly lists the underlying equations but then provides the some code without any context while also introducing additional variables (i.e. m1, m2). This ODE function needs some context, currently no equations or additional information are provided:

    function double_pendulum(xdot,x,p,t)
        xdot[1]=x[2]
        xdot[2]=-((g*(2*m₁+m₂)*sin(x[1])+m₂*(g*sin(x[1]-2*x[3])+2*(L₂*x[4]^2+L₁*x[2]^2*cos(x[1]-x[3]))*sin(x[1]-x[3])))/(2*L₁*(m₁+m₂-m₂*cos(x[1]-x[3])^2)))
        xdot[3]=x[4]
        xdot[4]=(((m₁+m₂)*(L₁*x[2]^2+g*cos(x[1]))+L₂*m₂*x[4]^2*cos(x[1]-x[3]))*sin(x[1]-x[3]))/(L₂*(m₁+m₂-m₂*cos(x[1]-x[3])^2))
    end
    
    

    The initial equations block only makes sense in the next subchapter (Poincaré section).

  • 14

    as for code readable, I have a little suggestion

    eg. lorenz

    function lorenz!(du,u,p,t)
        σ,ρ,β = p
      
        du[1] = σ*(u[2]-u[1])
        du[2] = u[1]*(ρ-u[3]) - u[2]
        du[3] = u[1]*u[2] - β*u[3]
    end
    

    modify to :

    function lorenz!(du,u,p,t)
        σ,ρ,β = p
        x,y,z = u
      
        du[1] = σ*(y-x)
        du[2] = x*(ρ-z) -y
        du[3] = x*y - β*z
    end
    

    better or worse?

  • 15

    SciMLTutorial does not install on Apple Silicon macOS system

    I am not able to install SciMLTutorial.jl on an Apple Silicon macOS system. It gets stuck in a loop.

    julia> using SciMLTutorials
    
    julia> SciMLTutorials.open_notebooks()
    [ Info: Weaving /Users/gong/.julia/packages/SciMLTutorials/3qOxd/src/../tutorials/DiffEqUncertainty/01-expectation_introduction.jmd
    ┌ Info: Instantiating
    └   folder = "DiffEqUncertainty"
      Activating project at `~/.julia/packages/SciMLTutorials/3qOxd/tutorials/DiffEqUncertainty`
    ┌ Warning: The active manifest file is an older format with no julia version entry. Dependencies may have been resolved with a different julia version.
    └ @ ~/.julia/packages/SciMLTutorials/3qOxd/tutorials/DiffEqUncertainty/Manifest.toml:0
       Installed DiffEqUncertainty ────────── v1.8.0
       Installed x265_jll ─────────────────── v3.0.0+3
       Installed AdvancedPS ───────────────── v0.2.2
       Installed MutableArithmetics ───────── v0.2.19
       Installed AbstractMCMC ─────────────── v3.2.0
       Installed libfdk_aac_jll ───────────── v0.1.6+4
       Installed GR_jll ───────────────────── v0.57.2+0
       Installed Quadrature ───────────────── v1.8.1
       Installed DifferentialEquations ────── v6.17.1
       Installed JpegTurbo_jll ────────────── v2.0.1+3
       Installed ConsoleProgressMonitor ───── v0.1.2
       Installed OffsetArrays ─────────────── v1.9.0
       Installed NonlinearSolve ───────────── v0.3.8
       Installed StatisticalTraits ────────── v1.1.0
       Installed LoggingExtras ────────────── v0.4.6
       Installed Preferences ──────────────── v1.2.2
       Installed StatsFuns ────────────────── v0.9.8
       Installed FFTW ─────────────────────── v1.4.1
       Installed NNlib ────────────────────── v0.7.19
       Installed NaturalSort ──────────────── v1.0.0
       Installed Opus_jll ─────────────────── v1.3.1+3
       Installed PDMats ───────────────────── v0.11.0
       Installed Unitful ──────────────────── v1.7.0
       Installed AxisArrays ───────────────── v0.4.3
       Installed TimerOutputs ─────────────── v0.5.9
       Installed LabelledArrays ───────────── v1.6.1
       Installed ZipFile ──────────────────── v0.9.3
       Installed StaticArrays ─────────────── v1.2.0
       Installed EarCut_jll ───────────────── v2.1.5+1
       Installed ProgressMeter ────────────── v1.6.2
       Installed Polyester ────────────────── v0.3.1
       Installed DiffEqPhysics ────────────── v3.9.0
       Installed PlotUtils ────────────────── v1.0.10
       Installed RecipesPipeline ──────────── v0.3.2
       Installed NearestNeighbors ─────────── v0.4.8
       Installed Sundials_jll ─────────────── v5.2.0+1
       Installed RecursiveArrayTools ──────── v2.11.4
       Installed InitialValues ────────────── v0.2.10
       Installed HTTP ─────────────────────── v0.9.9
       Installed InvertedIndices ──────────── v1.0.0
       Installed BFloat16s ────────────────── v0.1.0
       Installed GPUArrays ────────────────── v6.4.1
       Installed Cairo_jll ────────────────── v1.16.0+6
       Installed Fontconfig_jll ───────────── v2.13.1+14
       Installed RandomNumbers ────────────── v1.4.0
       Installed CheapThreads ─────────────── v0.2.5
       Installed IntelOpenMP_jll ──────────── v2018.0.3+2
       Installed Memoize ──────────────────── v0.4.4
       Installed ZygoteRules ──────────────── v0.2.1
       Installed Static ───────────────────── v0.2.4
       Installed EllipticalSliceSampling ──── v0.4.3
       Installed SLEEFPirates ─────────────── v0.6.20
       Installed Distances ────────────────── v0.10.3
       Installed LLVM ─────────────────────── v3.7.1
       Installed Missings ─────────────────── v1.0.0
       Installed FiniteDiff ───────────────── v2.8.0
       Installed EllipsisNotation ─────────── v1.1.0
       Installed SteadyStateDiffEq ────────── v1.6.2
       Installed FFMPEG ───────────────────── v0.4.0
       Installed Bijectors ────────────────── v0.9.2
       Installed IRTools ──────────────────── v0.4.2
       Installed MathOptInterface ─────────── v0.9.22
       Installed Functors ─────────────────── v0.2.1
       Installed ArnoldiMethod ────────────── v0.1.0
       Installed DiffEqBase ───────────────── v6.62.2
       Installed SciMLBase ────────────────── v1.13.4
       Installed JSON ─────────────────────── v0.21.1
       Installed BandedMatrices ───────────── v0.16.9
       Installed NLopt ────────────────────── v0.6.2
       Installed Qt5Base_jll ──────────────── v5.15.2+0
       Installed Hwloc_jll ────────────────── v2.4.1+0
       Installed Symbolics ────────────────── v0.1.25
       Installed Turing ───────────────────── v0.15.24
       Installed SpecialFunctions ─────────── v1.4.1
       Installed BenchmarkTools ───────────── v1.0.0
       Installed Bzip2_jll ────────────────── v1.0.6+5
       Installed ProgressLogging ──────────── v0.1.4
       Installed CodecBzip2 ───────────────── v0.7.2
       Installed SentinelArrays ───────────── v1.3.0
       Installed LibVPX_jll ───────────────── v1.9.0+1
       Installed IfElse ───────────────────── v0.1.0
       Installed DiffEqNoiseProcess ───────── v5.7.3
       Installed Zygote ───────────────────── v0.6.11
       Installed GlobalSensitivity ────────── v1.0.0
       Installed NaNMath ──────────────────── v0.3.5
       Installed OrdinaryDiffEq ───────────── v5.56.0
       Installed DefineSingletons ─────────── v0.1.1
       Installed MicroCollections ─────────── v0.1.0
       Installed ThreadingUtilities ───────── v0.4.4
       Installed DynamicPPL ───────────────── v0.10.20
       Installed DimensionalPlotRecipes ───── v1.2.0
       Installed Clustering ───────────────── v0.14.2
       Installed TranscodingStreams ───────── v0.9.5
       Installed FriBidi_jll ──────────────── v1.0.5+6
       Installed MonteCarloIntegration ────── v0.0.2
       Installed NLSolversBase ────────────── v7.8.0
       Installed SafeTestsets ─────────────── v0.0.1
       Installed MathProgBase ─────────────── v0.7.8
       Installed ZeroMQ_jll ───────────────── v4.3.2+6
       Installed GLFW_jll ─────────────────── v3.3.4+0
       Installed Random123 ────────────────── v1.3.1
       Installed Sobol ────────────────────── v1.5.0
       Installed x264_jll ─────────────────── v2020.7.14+2
       Installed Ratios ───────────────────── v0.4.0
       Installed LeftChildRightSiblingTrees ─ v0.1.2
       Installed FreeType2_jll ────────────── v2.10.1+5
       Installed MappedArrays ─────────────── v0.4.0
       Installed DataValues ───────────────── v0.4.13
       Installed ChainRules ───────────────── v0.7.65
       Installed DataStructures ───────────── v0.18.9
       Installed GR ───────────────────────── v0.57.4
       Installed StatsAPI ─────────────────── v1.0.0
       Installed Tracker ──────────────────── v0.2.16
       Installed ExprTools ────────────────── v0.1.3
       Installed Compat ───────────────────── v3.30.0
       Installed MCMCChains ───────────────── v4.12.0
       Installed BangBang ─────────────────── v0.3.30
       Installed SciMLTutorials ───────────── v0.9.0
       Installed ArrayLayouts ─────────────── v0.7.0
       Installed Libtask ──────────────────── v0.5.1
       Installed Scratch ──────────────────── v1.0.3
       Installed StatsPlots ───────────────── v0.14.21
       Installed StatsBase ────────────────── v0.33.8
       Installed ModelingToolkit ──────────── v5.16.0
       Installed DistributionsAD ──────────── v0.6.26
       Installed DiffRules ────────────────── v1.0.2
       Installed Cuba_jll ─────────────────── v4.2.1+0
       Installed Plots ────────────────────── v1.15.2
       Installed StochasticDiffEq ─────────── v6.34.1
       Installed Libtask_jll ──────────────── v0.4.3+0
       Installed ColorSchemes ─────────────── v3.12.1
       Installed ArgCheck ─────────────────── v2.1.0
       Installed NamedArrays ──────────────── v0.9.5
       Installed Conda ────────────────────── v1.5.2
       Installed Sundials ─────────────────── v4.4.3
       Installed QuasiMonteCarlo ──────────── v0.2.2
       Installed HCubature ────────────────── v1.5.0
       Installed CUDAKernels ──────────────── v0.1.0
       Installed Weave ────────────────────── v0.10.8
       Installed AbstractPPL ──────────────── v0.1.4
       Installed LatticeRules ─────────────── v0.0.1
       Installed Mustache ─────────────────── v1.0.10
       Installed AdvancedMH ───────────────── v0.6.0
       Installed AdvancedVI ───────────────── v0.1.3
       Installed DiffEqCallbacks ──────────── v2.16.1
       Installed Parsers ──────────────────── v1.1.0
       Installed JSONSchema ───────────────── v0.3.3
       Installed PrettyTables ─────────────── v1.0.1
       Installed Setfield ─────────────────── v0.7.0
       Installed Libtiff_jll ──────────────── v4.1.0+2
       Installed ParameterizedFunctions ───── v5.10.0
       Installed DelayDiffEq ──────────────── v5.31.0
       Installed MKL_jll ──────────────────── v2021.1.1+1
       Installed ScientificTypes ──────────── v1.1.2
       Installed Libffi_jll ───────────────── v3.2.2+0
       Installed ConstructionBase ─────────── v1.2.1
       Installed DiffEqFinancial ──────────── v2.4.0
       Installed DiffEqJump ───────────────── v6.14.2
       Installed Cuba ─────────────────────── v2.2.0
       Installed InplaceOps ───────────────── v0.3.0
       Installed Ogg_jll ──────────────────── v1.3.4+2
       Installed ChainRulesCore ───────────── v0.9.44
       Installed Combinatorics ────────────── v1.0.2
       Installed KernelDensity ────────────── v0.6.3
       Installed ArrayInterface ───────────── v3.1.15
       Installed Reexport ─────────────────── v1.0.0
       Installed TerminalLoggers ──────────── v0.1.3
       Installed LogExpFunctions ──────────── v0.2.4
       Installed MacroTools ───────────────── v0.5.6
       Installed KernelAbstractions ───────── v0.6.3
       Installed RangeArrays ──────────────── v0.3.2
       Installed FFTW_jll ─────────────────── v3.3.9+7
       Installed DiffEqSensitivity ────────── v6.45.0
       Installed CompositionsBase ─────────── v0.1.0
       Installed OpenSpecFun_jll ──────────── v0.5.4+0
       Installed MultivariateStats ────────── v0.8.0
       Installed YAML ─────────────────────── v0.4.6
       Installed SimpleTraits ─────────────── v0.9.3
       Installed Widgets ──────────────────── v0.6.3
       Installed DocStringExtensions ──────── v0.8.4
       Installed Wayland_protocols_jll ────── v1.18.0+4
       Installed RandomExtensions ─────────── v0.4.3
       Installed ExponentialUtilities ─────── v1.8.4
       Installed Transducers ──────────────── v0.4.65
       Installed NLopt_jll ────────────────── v2.7.0+0
       Installed Adapt ────────────────────── v3.3.0
       Installed AdvancedHMC ──────────────── v0.2.27
       Installed StableRNGs ───────────────── v1.0.0
       Installed LightGraphs ──────────────── v1.3.5
       Installed libass_jll ───────────────── v0.14.0+4
       Installed StrideArraysCore ─────────── v0.1.11
       Installed IntervalSets ─────────────── v0.5.3
       Installed Latexify ─────────────────── v0.15.5
       Installed VectorizationBase ────────── v0.20.11
       Installed Wayland_jll ──────────────── v1.17.0+4
       Installed MLJModelInterface ────────── v1.1.0
       Installed Arpack_jll ───────────────── v3.5.0+3
       Installed IterTools ────────────────── v1.3.0
       Installed RuntimeGeneratedFunctions ── v0.5.2
       Installed GPUCompiler ──────────────── v0.10.0
       Installed Cassette ─────────────────── v0.3.6
       Installed OpenSSL_jll ──────────────── v1.1.1+6
       Installed FFMPEG_jll ───────────────── v4.3.1+4
       Installed MultiScaleArrays ─────────── v1.8.1
       Installed GeometryBasics ───────────── v0.3.12
       Installed Tables ───────────────────── v1.4.2
       Installed VertexSafeGraphs ─────────── v0.1.2
       Installed Optim ────────────────────── v1.3.0
       Installed Interpolations ───────────── v0.13.2
       Installed Baselet ──────────────────── v0.1.1
       Installed DataAPI ──────────────────── v1.6.0
       Installed RecipesBase ──────────────── v1.1.1
       Installed TableOperations ──────────── v1.0.0
       Installed IterativeSolvers ─────────── v0.9.1
       Installed LAME_jll ─────────────────── v3.100.0+3
       Installed ReverseDiff ──────────────── v1.9.0
       Installed AbstractAlgebra ──────────── v0.16.0
       Installed Libiconv_jll ─────────────── v1.16.1+0
       Installed VersionParsing ───────────── v1.2.0
       Installed LatinHypercubeSampling ───── v1.8.0
       Installed LaTeXStrings ─────────────── v1.2.1
       Installed StructArrays ─────────────── v0.5.1
       Installed Arpack ───────────────────── v0.4.0
       Installed libvorbis_jll ────────────── v1.3.6+6
       Installed Glib_jll ─────────────────── v2.68.1+0
       Installed FillArrays ───────────────── v0.11.7
       Installed QuadGK ───────────────────── v2.4.1
       Installed ResettableStacks ─────────── v1.1.0
       Installed ForwardDiff ──────────────── v0.10.18
       Installed SymbolicUtils ────────────── v0.11.2
       Installed LoopVectorization ────────── v0.12.23
       Installed AxisAlgorithms ───────────── v1.0.0
       Installed Requires ─────────────────── v1.1.3
       Installed SplittablesBase ──────────── v0.1.13
       Installed FastBroadcast ────────────── v0.1.8
       Installed Distributions ────────────── v0.24.18
       Installed Trapz ────────────────────── v2.0.2
       Installed SortingAlgorithms ────────── v1.0.0
       Installed Parameters ───────────────── v0.12.2
       Installed RecursiveFactorization ───── v0.1.12
       Installed DiffEqGPU ────────────────── v1.12.0
       Installed WoodburyMatrices ─────────── v0.5.3
       Installed SparseDiffTools ──────────── v1.13.2
       Installed CUDA ─────────────────────── v2.6.3
      Downloaded artifact: Libffi
      Downloaded artifact: OpenSpecFun
      Downloaded artifact: Wayland_protocols
      Downloaded artifact: Libiconv
      Downloaded artifact: Glib
        Building Conda ────→ `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/299304989a5e6473d985212c28928899c74e9421/build.log`
        Building GR ───────→ `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/011458b83178ac913dc4eb73b229af45bdde5d83/build.log`
        Building Random123 → `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/7c6710c8198fd4444b5eb6a3840b7d47bd3593c5/build.log`
    Precompiling project...
      Progress [========================================>]  300/300
      ✓ VersionParsing
      ✓ Preferences
      ✓ Parsers
      ✓ Requires
      ✓ DocStringExtensions
      ✓ DataAPI
      ✓ YAML
      ✓ Scratch
      ✗ LLVM
      ✓ JLLWrappers
      ✓ Highlights
      ✓ Tables
      ✓ JSON
      ✓ libsodium_jll
      ✓ Libiconv_jll
      ✗ ChainRules
      ✓ VersionParsing
      ✓ Preferences
      ✓ Parsers
      ✓ Requires
      ✓ DocStringExtensions
      ✓ DataAPI
      ✓ YAML
      ✓ Scratch
      ✗ LLVM
      ✓ JLLWrappers
      ✓ Highlights
      ✓ Tables
      ✓ JSON
      ✓ libsodium_jll
      ✓ Libiconv_jll
      ✗ ChainRules
      ✗ GPUCompiler
      ✓ Conda
      ✓ Mustache
      ✓ ZeroMQ_jll
      ✗ FFTW
      ✗ CUDA
      ✓ Weave
      ✓ ZMQ
      ✗ IJulia
      ✗ CUDAKernels
      ✗ Zygote
      ✗ DistributionsAD
      ✗ KernelDensity
      ✗ GlobalSensitivity
      ✗ NLopt
      ✗ AdvancedVI
      ✗ DiffEqGPU
      ✗ Quadrature
      ✗ Turing
      ✗ SciMLTutorials
      ✗ StatsPlots
      ✗ DiffEqUncertainty
      ✗ DiffEqSensitivity
      280 dependencies successfully precompiled in 118 seconds (31 already precompiled)
      19 dependencies precompiled but different versions are currently loaded. Restart julia to access the new versions
      20 dependencies errored. To see a full report either run `import Pkg; Pkg.precompile()` or load the packages
        Building GR ───────→ `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/011458b83178ac913dc4eb73b229af45bdde5d83/build.log`
        Building Conda ────→ `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/299304989a5e6473d985212c28928899c74e9421/build.log`
        Building IJulia ───→ `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/d8b9c31196e1dd92181cd0f5760ca2d2ffb4ac0f/build.log`
        Building Random123 → `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/7c6710c8198fd4444b5eb6a3840b7d47bd3593c5/build.log`
    Precompiling project...
      4 dependencies successfully precompiled in 9 seconds (307 already precompiled, 20 skipped during auto due to previous errors)
    [ Info: Printing out `Pkg.status()`
          Status `~/.julia/packages/SciMLTutorials/3qOxd/tutorials/DiffEqUncertainty/Project.toml`
      [6e4b80f9] BenchmarkTools v1.0.0
      [8a292aeb] Cuba v2.2.0
      [071ae1c0] DiffEqGPU v1.12.0
      [41bf760c] DiffEqSensitivity v6.45.0
      [ef61062a] DiffEqUncertainty v1.8.0
      [0c46a032] DifferentialEquations v6.17.1
      [31c24e10] Distributions v0.24.18
      [f6369f11] ForwardDiff v0.10.18
      [5ab0869b] KernelDensity v0.6.3
      [c7f686f2] MCMCChains v4.12.0
      [76087f3c] NLopt v0.6.2
      [1dea7af3] OrdinaryDiffEq v5.56.0
      [91a5bcdd] Plots v1.15.2
      [67601950] Quadrature v1.8.1
      [30cb0354] SciMLTutorials v0.9.0
      [f3b207a7] StatsPlots v0.14.21
      [fce5fe82] Turing v0.15.24
    Building Notebook
    [ Info: Weaving /Users/gong/.julia/packages/SciMLTutorials/3qOxd/src/../tutorials/DiffEqUncertainty/02-AD_and_optimization.jmd
    ┌ Info: Instantiating
    └   folder = "DiffEqUncertainty"
      Activating project at `~/.julia/packages/SciMLTutorials/3qOxd/tutorials/DiffEqUncertainty`
    Precompiling project...
      4 dependencies successfully precompiled in 9 seconds (307 already precompiled, 20 skipped during auto due to previous errors)
    

    This is the part that loops again and again:

    [ Info: Printing out `Pkg.status()`
          Status `~/.julia/packages/SciMLTutorials/3qOxd/tutorials/advanced/Project.toml`
      [2169fc97] AlgebraicMultigrid v0.4.0
      [6e4b80f9] BenchmarkTools v1.0.0
      [052768ef] CUDA v2.6.3
      [2b5f629d] DiffEqBase v6.62.2
      [9fdde737] DiffEqOperators v4.26.0
      [0c46a032] DifferentialEquations v6.17.1
      [587475ba] Flux v0.12.1
      [961ee093] ModelingToolkit v5.17.3
      [2774e3e8] NLsolve v4.5.1
      [315f7962] NeuralPDE v3.10.1
      [1dea7af3] OrdinaryDiffEq v5.56.0
      [91a5bcdd] Plots v1.15.2
      [0bca4576] SciMLBase v1.13.4
      [30cb0354] SciMLTutorials v0.9.0
      [47a9eef4] SparseDiffTools v1.13.2
      [684fba80] SparsityDetection v0.3.4
      [789caeaf] StochasticDiffEq v6.34.1
      [c3572dad] Sundials v4.4.3
      [37e2e46d] LinearAlgebra
      [2f01184e] SparseArrays
    Building Notebook
    [ Info: Weaving /Users/gong/.julia/packages/SciMLTutorials/3qOxd/src/../tutorials/exercises/01-workshop_exercises.jmd
    ┌ Info: Instantiating
    └   folder = "exercises"
      Activating project at `~/.julia/packages/SciMLTutorials/3qOxd/tutorials/exercises`
       Installed RootedTrees ──── v1.0.0
       Installed DiffEqDevTools ─ v2.27.2
    Precompiling project...
      ✗ LLVM
      ✗ ChainRules
      ✗ FFTW
      ✗ SciMLTutorials
      ✗ KernelDensity
      ✗ GPUCompiler
      ✗ Zygote
      ✗ DistributionsAD
      ✗ CUDA
      ✗ GlobalSensitivity
      ✗ Flux
      ✗ CUDAKernels
      ✗ DiffEqSensitivity
      ✗ GalacticOptim
      ✗ DiffEqGPU
      ✗ DiffEqFlux
      7 dependencies successfully precompiled in 20 seconds (270 already precompiled)
      16 dependencies errored. To see a full report either run `import Pkg; Pkg.precompile()` or load the packages
        Building GR ───────→ `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/011458b83178ac913dc4eb73b229af45bdde5d83/build.log`
        Building Conda ────→ `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/299304989a5e6473d985212c28928899c74e9421/build.log`
        Building IJulia ───→ `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/d8b9c31196e1dd92181cd0f5760ca2d2ffb4ac0f/build.log`
        Building Random123 → `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/7c6710c8198fd4444b5eb6a3840b7d47bd3593c5/build.log`
        Building TimeZones → `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/960099aed321e05ac649c90d583d59c9309faee1/build.log`
    Precompiling project...
      4 dependencies successfully precompiled in 9 seconds (273 already precompiled, 16 skipped during auto due to previous errors)
    
    julia> versioninfo()
    Julia Version 1.7.1
    Commit ac5cc99908 (2021-12-22 19:35 UTC)
    Platform Info:
      OS: macOS (arm64-apple-darwin21.2.0)
      CPU: Apple M1 Pro
      WORD_SIZE: 64
      LIBM: libopenlibm
      LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
    
  • 16

    automatic differentiation not working in advanced examples?

    I'm trying to find an example of how to use automatic differentiation in more complex ODE rather than those vanilla examples in the tutorials. Simple ODEs are quite easy to figure out, but I always have the trouble with the Dual Number error in complex ODE models. Seems my ODE functions are always not generic enough. And I'd like to learn how to make it work.

    I was trying to modify the Beeler-Reuter Model (the CPU version) in the advanced tutorial, as the original code doesn't work with autodiff=true since it only accepts Float. I finally manage to make it work, but I have many problems that I don't understand. Here is the modified working code:

    using DifferentialEquations, Sundials, ForwardDiff
    
    const v0 = -84.624
    const v1 = 10.0
    const C_K1 = 1.0f0
    const C_x1 = 1.0f0
    const C_Na = 1.0f0
    const C_s = 1.0f0
    const D_Ca = 0.0f0
    const D_Na = 0.0f0
    const g_s = 0.09f0
    const g_Na = 4.0f0
    const g_NaC = 0.005f0
    const ENa = 50.0f0 + D_Na
    const γ = 0.5f0
    const C_m = 1.0f0
    
    # make this parametric
    mutable struct BeelerReuterCpu{T} <: Function
        t::T              # the last timestep time to calculate Δt
        diff_coef::T      # the diffusion-coefficient (coupling strength)
    
        C::Array{T,2}   # intracellular calcium concentration
        M::Array{T,2}     # sodium current activation gate (m)
        H::Array{T,2}    # sodium current inactivation gate (h)
        J::Array{T,2}     # sodium current slow inactivaiton gate (j)
        D::Array{T,2}     # calcium current activaiton gate (d)
        F::Array{T,2}     # calcium current inactivation gate (f)
        XI::Array{T,2}    # inward-rectifying potassium current (iK1)
        Δu::Array{T,2}    # place-holder for the Laplacian
    
        function BeelerReuterCpu(u0::Array{T,2}, diff_coef::T) where {T}
            self = new{T}()
    
            ny, nx = size(u0)
            self.t = 0.0
            self.diff_coef = diff_coef
    
            self.C = fill(0.0001f0, (ny, nx))
            self.M = fill(0.01f0, (ny, nx))
            self.H = fill(0.988f0, (ny, nx))
            self.J = fill(0.975f0, (ny, nx))
            self.D = fill(0.003f0, (ny, nx))
            self.F = fill(0.994f0, (ny, nx))
            self.XI = fill(0.0001f0, (ny, nx))
            self.Δu = zeros(ny, nx)
    
            return self
        end
    end
    
    
    # 5-point stencil
    function laplacian(Δu, u) 
    
        n1, n2 = size(u)
    
        # internal nodes
        for j = 2:n2-1
            for i = 2:n1-1
                @inbounds Δu[i, j] = u[i+1, j] + u[i-1, j] + u[i, j+1] + u[i, j-1] - 4 * u[i, j]
            end
        end
    
        # left/right edges
        for i = 2:n1-1
            @inbounds Δu[i, 1] = u[i+1, 1] + u[i-1, 1] + 2 * u[i, 2] - 4 * u[i, 1]
            @inbounds Δu[i, n2] = u[i+1, n2] + u[i-1, n2] + 2 * u[i, n2-1] - 4 * u[i, n2]
        end
    
        # top/bottom edges
        for j = 2:n2-1
            @inbounds Δu[1, j] = u[1, j+1] + u[1, j-1] + 2 * u[2, j] - 4 * u[1, j]
            @inbounds Δu[n1, j] = u[n1, j+1] + u[n1, j-1] + 2 * u[n1-1, j] - 4 * u[n1, j]
        end
    
        # corners
        @inbounds Δu[1, 1] = 2 * (u[2, 1] + u[1, 2]) - 4 * u[1, 1]
        @inbounds Δu[n1, 1] = 2 * (u[n1-1, 1] + u[n1, 2]) - 4 * u[n1, 1]
        @inbounds Δu[1, n2] = 2 * (u[2, n2] + u[1, n2-1]) - 4 * u[1, n2]
        @inbounds Δu[n1, n2] = 2 * (u[n1-1, n2] + u[n1, n2-1]) - 4 * u[n1, n2]
    end
    
    @inline function rush_larsen(g, α, β, Δt)
        inf = α / (α + β)
        τ = 1.0f0 / (α + β)
        return clamp(g + (g - inf) * expm1(-Δt / τ), 0.0f0, 1.0f0)
    end
    
    function update_M_cpu(g, v, Δt)
        # the condition is needed here to prevent NaN when v == 47.0
        α = isapprox(v, 47.0f0) ? 10.0f0 : -(v + 47.0f0) / (exp(-0.1f0 * (v + 47.0f0)) - 1.0f0)
        β = (40.0f0 * exp(-0.056f0 * (v + 72.0f0)))
        return rush_larsen(g, α, β, Δt)
    end
    
    function update_H_cpu(g, v, Δt)
        α = 0.126f0 * exp(-0.25f0 * (v + 77.0f0))
        β = 1.7f0 / (exp(-0.082f0 * (v + 22.5f0)) + 1.0f0)
        return rush_larsen(g, α, β, Δt)
    end
    
    function update_J_cpu(g, v, Δt)
        α = (0.55f0 * exp(-0.25f0 * (v + 78.0f0))) / (exp(-0.2f0 * (v + 78.0f0)) + 1.0f0)
        β = 0.3f0 / (exp(-0.1f0 * (v + 32.0f0)) + 1.0f0)
        return rush_larsen(g, α, β, Δt)
    end
    
    function update_D_cpu(g, v, Δt)
        α = γ * (0.095f0 * exp(-0.01f0 * (v - 5.0f0))) / (exp(-0.072f0 * (v - 5.0f0)) + 1.0f0)
        β = γ * (0.07f0 * exp(-0.017f0 * (v + 44.0f0))) / (exp(0.05f0 * (v + 44.0f0)) + 1.0f0)
        return rush_larsen(g, α, β, Δt)
    end
    
    function update_F_cpu(g, v, Δt)
        α = γ * (0.012f0 * exp(-0.008f0 * (v + 28.0f0))) / (exp(0.15f0 * (v + 28.0f0)) + 1.0f0)
        β = γ * (0.0065f0 * exp(-0.02f0 * (v + 30.0f0))) / (exp(-0.2f0 * (v + 30.0f0)) + 1.0f0)
        return rush_larsen(g, α, β, Δt)
    end
    
    function update_XI_cpu(g, v, Δt)
        α = (0.0005f0 * exp(0.083f0 * (v + 50.0f0))) / (exp(0.057f0 * (v + 50.0f0)) + 1.0f0)
        β = (0.0013f0 * exp(-0.06f0 * (v + 20.0f0))) / (exp(-0.04f0 * (v + 20.0f0)) + 1.0f0)
        return rush_larsen(g, α, β, Δt)
    end
    
    function update_C_cpu(g, d, f, v, Δt)
        ECa = D_Ca - 82.3f0 - 13.0278f0 * log(g)
        kCa = C_s * g_s * d * f
        iCa = kCa * (v - ECa)
        inf = 1.0f-7 * (0.07f0 - g)
        τ = 1.0f0 / 0.07f0
        return g + (g - inf) * expm1(-Δt / τ)
    end
    
    function update_gates_cpu(u, XI, M, H, J, D, F, C, Δt)
        # let Δt = Float32(Δt) # remove the Let 
            n1, n2 = size(u)
            for j = 1:n2
                for i = 1:n1
                    v = u[i, j]
    
                    XI[i, j] = update_XI_cpu(XI[i, j], v, Δt)
                    M[i, j] = update_M_cpu(M[i, j], v, Δt)
                    H[i, j] = update_H_cpu(H[i, j], v, Δt)
                    J[i, j] = update_J_cpu(J[i, j], v, Δt)
                    D[i, j] = update_D_cpu(D[i, j], v, Δt)
                    F[i, j] = update_F_cpu(F[i, j], v, Δt)
    
                    C[i, j] = update_C_cpu(C[i, j], D[i, j], F[i, j], v, Δt)
                end
            end
        # end
    end
    
    # iK1 is the inward-rectifying potassium current
    function calc_iK1(v)
        ea = exp(0.04f0 * (v + 85.0f0))
        eb = exp(0.08f0 * (v + 53.0f0))
        ec = exp(0.04f0 * (v + 53.0f0))
        ed = exp(-0.04f0 * (v + 23.0f0))
        return 0.35f0 * (
            4.0f0 * (ea - 1.0f0) / (eb + ec) +
            0.2f0 * (isapprox(v, -23.0f0) ? 25.0f0 : (v + 23.0f0) / (1.0f0 - ed))
        )
    end
    
    # ix1 is the time-independent background potassium current
    function calc_ix1(v, xi)
        ea = exp(0.04f0 * (v + 77.0f0))
        eb = exp(0.04f0 * (v + 35.0f0))
        return xi * 0.8f0 * (ea - 1.0f0) / eb
    end
    
    # iNa is the sodium current (similar to the classic Hodgkin-Huxley model)
    function calc_iNa(v, m, h, j)
        return C_Na * (g_Na * m^3 * h * j + g_NaC) * (v - ENa)
    end
    
    # iCa is the calcium current
    function calc_iCa(v, d, f, c)
        ECa = D_Ca - 82.3f0 - 13.0278f0 * log(c)    # ECa is the calcium reversal potential
        return C_s * g_s * d * f * (v - ECa)
    end
    
    function update_du_cpu(du, u, XI, M, H, J, D, F, C)
        n1, n2 = size(u)
    
        for j = 1:n2
            for i = 1:n1
                v = u[i, j]
    
                # calculating individual currents
                iK1 = calc_iK1(v)
                ix1 = calc_ix1(v, XI[i, j])
                iNa = calc_iNa(v, M[i, j], H[i, j], J[i, j])
                iCa = calc_iCa(v, D[i, j], F[i, j], C[i, j])
    
                # # total current
                I_sum = iK1 + ix1 + iNa + iCa
    
                # # the reaction part of the reaction-diffusion equation
                du[i, j] = -I_sum / C_m
            end
        end
    end
    
    function (f::BeelerReuterCpu)(du, u, p, t)
        Δt = t - f.t
    
        if Δt != 0 || t == 0
            update_gates_cpu(
                u,
                eltype(u).(f.XI), # type conversion
                eltype(u).(f.M),
                eltype(u).(f.H),
                eltype(u).(f.J),
                eltype(u).(f.D),
                eltype(u).(f.F),
                eltype(u).(f.C),
                eltype(u)(Δt),
            )
            f.t = t
        end
    
        laplacian(eltype(u).(f.Δu), u) # type conversion
    
        # calculate the reaction portion
        update_du_cpu(du, u, f.XI, f.M, f.H, f.J, f.D, f.F, f.C)
    
        # ...add the diffusion portion
        du .+= f.diff_coef .* f.Δu
        du
    end
    
    N = 10;
    u0 = fill(v0, (N, N));
    u0[5:6,5:6] .= v1;   # a small square in the middle of the domain
    
    using Plots
    heatmap(u0)
    
    deriv_cpu = BeelerReuterCpu(u0, 1.0);
    prob = ODEProblem(deriv_cpu, u0, (0.0, 5.0))
    # Lower tolerances to show the methods converge to the same value
    @time sol=solve(prob, TRBDF2(autodiff = true), saveat = 100.0)
    

    The main changes are :

    1. Make the struct BeelerReuterCpu parametric.
    2. In the (f::BeelerReuterCpu)(du, u, p, t) function, when calling update_gates_cpu and laplacian functions, the aruguments are converted to the type of u , like eltype(u).(f.Δu).

    I made these changes otherwise I get the Dual Number error

    ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, BeelerReuterCpu{Float64}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, SciMLBase.NullParameters}, Float64}, Float64, 12})
    Closest candidates are:
      (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
      (::Type{T})(::T) where T<:Number at boot.jl:760
      (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
    

    My questions are:

    1. why type conversion is needed when calling update_gates_cpu and laplacian but not update_du_cpu inside (f::BeelerReuterCpu)(du, u, p, t)?
    2. I used parametric type in BeelerReuterCpu to make the type consistent between u and the fields in BeelerReuterCpu . I assume that's necessary for the Dual Number to work, but obviously that's not enough. I suppose that's because the struct is not created dynamically? Is there a way to do that.
    3. In the Documentation it is said DiffEqBase.dualcache is needed to avoid the Dual Number error when using cache. The cache example given there was very straightforward. But I wonder what exactly is considered as cache. Does it include any intermediary calculation variables?

    I hope there could be an example of complex ODE with automatic differentiation. All the complex ODE examples in the tutorials and benchmarks only accept Float and therefore doesn't work for automatic differentiation. Thanks!