make nonllinchem an ODEProblem with analytic

Author: ChrisRackauckas

src/ode/filament_prob.jl

@@ -158,7 +158,7 @@ function (f::FilamentCache)(dr, r, p, t)
     mul!(S1, A, r)
     S1 .+= F
     mul!(S2, P, S1)
-    copy!(dr, S2)
+    copyto!(dr, S2)
     return dr
 end
 
@@ -183,10 +183,10 @@ end
 function projection!(f::FilamentCache)
     # implement P[:] = I - J'/(J*J')*J in an optimized way to avoid temporaries
     J, P, J_JT, J_JT_LDLT, P0 = f.P.J, f.P.P, f.P.J_JT, f.P.J_JT_LDLT, f.P.P0
-    A_mul_Bt!(J_JT, J, J)
+    mul!(J_JT, J, J')
     LDLt_inplace!(J_JT_LDLT, J_JT)
-    A_ldiv_B!(P0, J_JT_LDLT, J)
-    At_mul_B!(P, P0, J)
+    ldiv!(P0, J_JT_LDLT, J)
+    mul!(P', P0, J)
     subtract_from_identity!(P)
     nothing
 end

…of_reactions_with_an_analytical_solution src/ode/nonlinchem.jlsrc/ode/ODE_Nonlinear_system_of_reactions_with_an_analytical_solution renamed to src/ode/nonlinchem.jl src/ode/ODE_Nonlinear_system_of_reactions_with_an_analytical_solution renamed to src/ode/nonlinchem.jl

@@ -1,28 +1,12 @@
-using DifferentialEquations
-using Plots
-using SpecialFunctions
-
 function nonLinChem(dy,y,p,t)
  dy[1] = -y[1]
  dy[2] = y[1]-(y[2])^2
  dy[3] = (y[2])^2
 end
 y0 = [1.0;0.0;0.0]
 tspan = (0.0,20)
-prob = ODEProblem(nonLinChem,y0,tspan)
-sol = solve(prob)
-
-Analytical(t) = [exp(-t);
+nlc_analytic(u0,p,t) = [exp(-t);
     (2sqrt(exp(-t))besselk(1,2sqrt(exp(-t)))-2besselk(1,2)/besseli(1,2)*sqrt(exp(-t))besseli(1,2sqrt(exp(-t))))/(2besselk(0,2sqrt(exp(-t)))+(2besselk(1,2)/besseli(1,2))besseli(0,2sqrt(exp(-t))));
     -exp(-t)+1+(-2sqrt(exp(-t))*besselk(1,2sqrt(exp(-t)))+sqrt(exp(-t))*besseli(1,2sqrt(exp(-t)))*2besselk(1,2)/besseli(1,2))/(2besselk(0,2sqrt(exp(-t)))+2besselk(1,2)/besseli(1,2)*besseli(0,2sqrt(exp(-t))))]
-
-t=0:.1:20
-p = 1
-plot(sol,vars=(0,p))
-plot!(t,(x->Analytical(x)[p]).(t))
-p = 2
-plot!(sol,vars=(0,p))
-plot!(t,(x->Analytical(x)[p]).(t))
-p = 3
-plot!(sol,vars=(0,p))
-plot!(t,(x->Analytical(x)[p]).(t))
+nonLinChem_f = ODEFunction(nonLinChem,analytic = nlc_analytic)
+prob_ode_nonlinchem = ODEProblem(nonLinChem,y0,tspan)

src/ode/ode_premade_problems.jl

@@ -10,10 +10,12 @@ export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear,
        prob_ode_lorenz, prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades,
        prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego,
        prob_ode_hires, prob_ode_pollution, prob_ode_filament,
+       prob_ode_nonlinchem,
        SolverDiffEq, filament_prob
 
 include("ode_linear_prob.jl")
 include("ode_simple_nonlinear_prob.jl")
 include("brusselator_prob.jl")
 include("pollution_prob.jl")
 include("filament_prob.jl")
+include("nonlinchem.jl")