diff --git a/Project.toml b/Project.toml index fac7af1..259a23b 100644 --- a/Project.toml +++ b/Project.toml @@ -4,20 +4,30 @@ version = "0.1.0" authors = ["Olivier Cots "] [deps] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" +DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c" +ODEInterfaceDiffEq = "09606e27-ecf5-54fc-bb29-004bd9f985bf" +ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] +DiffEqDevTools = "2.48.0" DifferentiationInterface = "0.7.9" +DualNumbers = "0.6.9" ForwardDiff = "1.2.2" LaTeXStrings = "1.4.0" LinearAlgebra = "1.12.0" +ODEInterface = "0.5.0" +ODEInterfaceDiffEq = "3.13.5" +ODEProblemLibrary = "1.2.2" OrdinaryDiffEq = "6.102.1" Plots = "1.41.1" SciMLSensitivity = "7.90.0" diff --git a/article/figures/chris.jl b/article/figures/chris.jl new file mode 100644 index 0000000..1b062db --- /dev/null +++ b/article/figures/chris.jl @@ -0,0 +1,56 @@ +############################################################################### +# From https://github.com/SciML/ODE.jl/blob/master/src/ODE.jl +# estimator for initial step based on book +# "Solving Ordinary Differential Equations I" by Hairer et al., p.169 +# +#function hinit(F, x0, t0::T, tend, p, reltol, abstol) where T +function hinit(F, x0, t0, tend, p, reltol, abstol) + # Returns first step, direction of integration and F evaluated at t0 + tdir = sign(tend-t0) + tdir==0 && error("Zero time span") + tau = max(reltol*norm(x0, Inf), abstol) + d0 = norm(x0, Inf)/tau + f0 = F(t0, x0) + d1 = norm(f0, Inf)/tau + if d0 < 1e-5 || d1 < 1e-5 + h0 = 1e-6 + else + h0 = 0.01*(d0/d1) + end + # perform Euler step + x1 = x0 + tdir*h0*f0 + f1 = F(t0 + tdir*h0, x1) + # estimate second derivative + d2 = norm(f1 - f0, Inf)/(tau*h0) + if max(d1, d2) <= 1e-15 + h1 = max((10)^(-6), (10)^(-3)*h0) + else + pow = -(2 + log10(max(d1, d2)))/(p + 1) + h1 = 10^pow + end + + return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0 +end + + +A(λ) = [λ[1] 0 0 ; 0 λ[2] 0 ; 0 0 λ[1]-λ[2]] +fun1(t,x) = A(λ)*x + +fun2(t,xX) = [fun1(xX[:,1],λ,t) fun1(xX[:,2:end],λ,t)+[xX[1,1] 0 ; 0 xX[2,1] ; xX[3,1] xX[3,1]] ] + +t0 = 0. ; tend = 2.; +p = 2 +x0 = [0.,1.,1] +reltol = 1.e-4; abstol = 1.e-4 + +hinit(fun1, x0, t0, tend, p, reltol, abstol) + +xX0 = [funx0(λ) [0 1 ; 0 0 ; 0 0]] +hinit(fun2, xX0, t0, tend, p, reltol, abstol) +RelTol = reltol*ones(n,p+1) +AbsTol = RelTol + +hinit(fun2, xX0, t0, tend, p, RelTol, AbsTol) + +RelTol = [reltol*ones(n,1) Inf*ones(n,p)]/sqrt(p+1) +AbsTol = [abstol*ones(n,1) Inf*ones(n,p)]/sqrt(p+1) diff --git a/article/figures/test-dopri-DP5-julia.jl b/article/figures/test-dopri-DP5-julia.jl new file mode 100644 index 0000000..181364c --- /dev/null +++ b/article/figures/test-dopri-DP5-julia.jl @@ -0,0 +1,48 @@ +using OrdinaryDiffEq, Test, DiffEqDevTools, + ODEInterface, ODEInterfaceDiffEq + +import ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear + +prob = prob_ode_linear +#= +prob.tspan +prob.f +prob.u0 +prob.p +prob.kwargs +prob.problem_type +=# + +sol = solve(prob, DP5()) + +sol2 = solve(prob, dopri5()) +println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2]) +@test sol.t[2] ≈ sol2.t[2] + +prob = prob_ode_2Dlinear +sol = solve(prob, DP5(), internalnorm = (u, t) -> sqrt(sum(abs2, u))) + +# Change the norm due to error in dopri5.f +sol2 = solve(prob, dopri5()) +println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2]) +@test sol.t[2] ≈ sol2.t[2] + +prob = deepcopy(prob_ode_linear) +prob +#prob2 = ODEProblem(prob.f, prob.u0, (1.0, 0.0), 1.01) +prob2 = ODEProblem(prob.f, prob.u0, (0.0, 1.0), 1.01) +sol = solve(prob2, DP5()) + +sol2 = solve(prob2, dopri5()) +println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2]) +@test sol.t[2] ≈ sol2.t[2] + +prob = deepcopy(prob_ode_2Dlinear) +#prob2 = ODEProblem(prob.f, prob.u0, (1.0, 0.0), 1.01) +prob2 = ODEProblem(prob.f, prob.u0, (0.0, 1.0), 1.01) +sol = solve(prob2, DP5(), internalnorm = (u, t) -> norm(u)/sqrt(length(u)))#sqrt(sum(abs2, u))) + +# Change the norm due to error in dopri5.f +sol2 = solve(prob2, dopri5()) +println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2]) +@test sol.t[2] ≈ sol2.t[2] \ No newline at end of file diff --git a/article/figures/test_myode43.jl b/article/figures/test_myode43.jl new file mode 100644 index 0000000..5dece0a --- /dev/null +++ b/article/figures/test_myode43.jl @@ -0,0 +1,56 @@ +include("../../src/myode43/myode43.jl") + + +# ẋ = [λ₁x₁ ; λ₂x₂ ; (λ₁-λ₂)x₃ +# x(t0) = x₀(λ) = [λ₂, 1, 1] +# x(t,t0,x_0(λ),λ) = diag(exp(λ₁t),exp(λ₂t),exp((λ₁-λ₂)t))x₀(λ) +# (∂x(t,t0,x_0(λ),λ)/∂x0) = diag(exp(λ₁t),exp(λ₂t),exp((λ₁-λ₂)t)) +# ∂x(t,t0,x_0(λ),λ)/∂λ = [λ₂texp(λ₁t) exp(λ₁t) ; 0 texp(λ₂t) ; texp((λ₁-λ₂)t)) -texp((λ₁-λ₂)t))] +# +# ∂x0_flow + +A(λ) = [λ[1] 0 0 ; 0 λ[2] 0 ; 0 0 λ[1]-λ[2]] +fun1(x,λ,t) = A(λ)*x + +fun2(xX,λ,t) = [fun1(xX[:,1],λ,t) fun1(xX[:,2:end],λ,t)+[xX[1,1] 0 ; 0 xX[2,1] ; xX[3,1] xX[3,1]] ] + + + +t0 = 0. ; tf = 2.; tspan = (t0,tf); +λ = [1.,2]; p = length(λ); +funx0(λ) = [λ[2],1.,1]; +x0 = funx0(λ) +n = length(x0) +reltol = 1.e-4; abstol = 1.e-4 +myode43(fun1,x0,λ,tspan,reltol,abstol) + +xX0 = [x0 [0 1 ; 0 0 ; 0 0]] + +myode43(fun2,xX0,λ,tspan,reltol,abstol) + +RelTol = reltol*ones(n,p+1) +AbsTol = RelTol + +myode43(fun2,xX0,λ,tspan,RelTol,AbsTol) + +RelTol = [reltol*ones(n,1) Inf*ones(n,p)]/sqrt(p+1) +AbsTol = [abstol*ones(n,1) Inf*ones(n,p)]/sqrt(p+1) + +inith(fun1,x0,λ,t0,abstol,reltol) + +inith(fun2,xX0,λ,t0,abstol,reltol) +RelTol = reltol*ones(n,p+1) +AbsTol = RelTol +inith(fun2,xX0,λ,t0,AbsTol,RelTol) + + +epsi = eps() +epsi = 0 +xX0 = [x0 [epsi 1 ; epsi epsi ; epsi epsi]] +my_Inf = prevfloat(typemax(Float64)) +#my_Inf = Inf +RelTol = [reltol*ones(n,1) my_Inf*ones(n,p)]/sqrt(p+1) +AbsTol = [abstol*ones(n,1) my_Inf*ones(n,p)]/sqrt(p+1) +inith(fun2,xX0,λ,t0,AbsTol,RelTol) +myode43(fun2,xX0,λ,tspan,RelTol,AbsTol) + diff --git a/article/figures/test_zygote.jl b/article/figures/test_zygote.jl new file mode 100644 index 0000000..d9080a3 --- /dev/null +++ b/article/figures/test_zygote.jl @@ -0,0 +1,93 @@ +using OrdinaryDiffEq +using LinearAlgebra +using Plots +using LaTeXStrings +using ForwardDiff: ForwardDiff +using Zygote: Zygote +#using Mooncake # plante à l'instalation +using SciMLSensitivity +using DifferentiationInterface +using DataFrames + +include("../../src/CTDiffFlow.jl") +using .CTDiffFlow +include("../../src/myode43/myode43.jl") + +function my_euler(fun,t0,x0,tf,λ,N) + ti = t0; xi = x0 + h = (tf-t0)/N + for i in 1:N + xi += h*fun(xi,λ,ti) + ti = ti+h + end + return xi +end + + +# Definition of the edo +A(λ) = [λ[1] 0 0 ; 0 λ[2] 0 ; 0 0 λ[1]-λ[2]] +fun1(x,λ,t) = A(λ)*x +funx0(λ) = [λ[2],1.,1]; +N = 10 +t0 = 0. ; tf = 1.; tspan = (t0,tf); +λ = [1.,2]; +reltol = 1.e-8 +abstol = reltol +# ∂λ_flow(λ) +sol_∂λ_flow = [λ[2]*exp(λ[1]*tf) exp(λ[1]*tf) + 0. tf*exp(λ[2]*tf) + tf*exp((λ[1]-λ[2])*tf) -tf*exp((λ[1]-λ[2])*tf)] + + +#dict_Zygote = Dict(:∂λ_sol => sol_∂λ_flow) +#dict_ForwardDiff = Dict(:∂λ_sol => sol_∂λ_flow) + +df = DataFrame(AutoDiff = String[], algo=String[], Jacobian = Matrix{<:Real}[]) +push!(df, ("solution", "", sol_∂λ_flow)) + +# Zygote with my_euler +_flow(λ) = my_euler(fun1,t0,funx0(λ),tf,λ,N) +#dict_Zygote[:my_euler] = jacobian(_flow,AutoZygote(),λ) +#dict_ForwardDiff[:my_euler] = jacobian(_flow,AutoForwardDiff(),λ) +push!(df, ("Zygote", "my-euler", jacobian(_flow,AutoZygote(),λ))) +push!(df, ("ForwardDiff", "my-euler", jacobian(_flow,AutoForwardDiff(),λ))) +# Zygote with numerical integration +function _flow_int(λ) + ivp = ODEProblem(fun1, funx0(λ), (t0,tf), λ) + #algo = get(ode_kwargs, :alg, Tsit5()) + #println("algo = ", algo) + + sol = solve(ivp, alg = Euler(),reltol=reltol,abstol=abstol,adaptive=false, dt = (tf-t0)/N) + return sol.u[end] +end + +#dict_Zygote[:Euler] = jacobian(_flow_int,AutoZygote(),λ) +#dict_ForwardDiff[:Euler] = jacobian(_flow_int,AutoForwardDiff(),λ) +push!(df, ("Zygote", "Euler", jacobian(_flow_int,AutoZygote(),λ))) +push!(df, ("ForwardDiff", "Euler", jacobian(_flow_int,AutoForwardDiff(),λ))) +# Zygote with ∂λ_flow +∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1,t0,funx0,tf,λ; backend=AutoZygote()) +#dict_Zygote[:CTDiffFlowZygoteEuler] = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N)#,print_times=false) +push!(df, ("Zygote", "CTDiffFlow-Euler", ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N))) + +# ForwardDiff with ∂λ_flow +∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1,t0,funx0,tf,λ; backend=AutoForwardDiff()) +#dict_ForwardDiff[:CTDiffFlowForwardDiffEuler] = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N)#,print_times=false) +push!(df, ("ForwardDiff", "my-euler", ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N))) + +function _flow(λ) + T, X = myode43(fun1,funx0(λ),λ,(t0,tf),abstol,reltol) + return X[end] +end + +#dict_Zygote[:myode43Zygote] = jacobian(_flow,AutoZygote(),λ) +#dict_ForwardDiff[:myode43ForwardDiff] = jacobian(_flow,AutoForwardDiff(),λ) +push!(df, ("Zygote", "my-ode43", [NaN;;])) +push!(df, ("ForwardDiff", "my-ode43", jacobian(_flow,AutoForwardDiff(),λ))) + +jacobian(_flow,AutoForwardDiff(),λ) + +latexify(df,arraystyle=:pmatrix,env=:tabular) + +sol_Zygote = df[in.(df.AutoDiff, Ref(["solution","Zygote"])),:] +sol_ForwardDiff = df[in.(df.AutoDiff, Ref(["solution", "ForwardDiff"])),:] diff --git a/article/figures/time-steps-bruss.jl b/article/figures/time-steps-bruss.jl new file mode 100644 index 0000000..f9b7684 --- /dev/null +++ b/article/figures/time-steps-bruss.jl @@ -0,0 +1,102 @@ + +using Pkg +Pkg.activate(".") +using OrdinaryDiffEq +using LinearAlgebra +using Plots +using LaTeXStrings +using ForwardDiff: ForwardDiff +using Zygote: Zygote +#using Mooncake # plante à l'instalation +using SciMLSensitivity +using DifferentiationInterface + +include("../../src/CTDiffFlow.jl") +using .CTDiffFlow +# +# Definition of the second member +function bruss(x,par,t) + λ = par[1] + x₁ = x[1] + x₂ = x[2] + return [1+x₁^2*x₂-(λ+1)*x₁ , λ*x₁-x₁^2*x₂] +end + +t0 = 0.; tf = 1. +tspan = (t0,tf) +λ = [3.] +p = length(λ) +x01 = 1.3 +funx0(λ) = [x01, λ[1]] +tol = 1.e-4 +n = 2 + +plt1 = plot(); plt2 = plot() +algo = Tsit5() +reltol = 1.e-4; abstol = 1.e-4 + +# Flow +ivp = ODEProblem(bruss, funx0(λ), (t0,tf), λ) +sol_ivp = solve(ivp, alg=algo; reltol=reltol,abstol=abstol) +println("Flow") +println("Times for the initial flow T = ") +println(sol_ivp.t) +dict_T = Dict(:ivp => sol_ivp.t) +dt = sol_ivp.t[2] + +# ------------------------------- +println("Vatiational equation") +# ------------------------------ + +∂λ_flow_var = CTDiffFlow.build_∂λ_flow_var(bruss,t0,funx0,tf,λ) + +#println(∂λ_flow_var(t0,funx0,tf,λ;reltol=reltol,abstol=abstol)) +sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=reltol,abstol=abstol) +sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=reltol,abstol=abstol,internalnorm = (u,t)->norm(u)) +println("Times default values= ", sol_var.t) +dict_T[:var_reltol] = sol_var.t +my_Inf = prevfloat(typemax(Float64)) +RelTol = [reltol*ones(n,1) Inf*ones(n,1)]/sqrt(p+1) +AbsTol = [abstol*ones(n,1) Inf*ones(n,1)]/sqrt(p+1) +sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=RelTol,abstol=AbsTol) +sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=RelTol,abstol=AbsTol, dt=dt)#,internalnorm = (u,t)->norm(u) ) +println("RelTol = ", RelTol) +println("Times = ", sol_var.t) +println(sol_ivp.t - sol_var.t) +dict_T[:var_RelTol] = sol_var.t + +# ne fonctionne pas +#@btime ∂λ_flow_var(t0,funx0,tf,λ;reltol=reltol,abstol=abstol) + +println("Automatic differentiation") +println("with ForwardDiff") +∂λ_flow = CTDiffFlow.build_∂λ_flow(bruss,t0,funx0,tf,λ) +sol_diff_auto_flow, T = ∂λ_flow(t0,funx0,tf,λ;reltol=RelTol,abstol=AbsTol,print_times=true) +println("Times : ", T) +dict_T[:diff_auto_ForwardDiff] = T +#sol_diff_auto_flow, T = ∂λ_flow(tspan,funx0,λ;reltol=RelTol,abstol=AbsTol,print_times=true) +#= +println("internalnorm = (u,t)->norm(u)") +println(∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,internalnorm = (u,t)->norm(u),print_times=true)) +#@btime ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol) +=# +println("With Zygote") +# with Zygote +∂λ_flow = CTDiffFlow.build_∂λ_flow(bruss,t0,funx0,tf,λ; backend=AutoZygote()) +sol_diff_auto_flow_tf, T = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,print_times=true) +dict_T[:diff_auto_Zygote] = T +println(sol_diff_auto_flow_tf-sol_var) +#sol_diff_auto_flow,T = ∂λ_flow(tspan,funx0,λ;reltol=reltol,abstol=abstol,print_times=true) +#println(sol_diff_auto_flow) +#reshape(sol_diff_auto_flow,2,7) + +#@btime ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol) + +println("With Mooncake") +# plante +#∂λ_flow = CTDiffFlow.build_∂λ_flow(bruss,t0,funx0,tf,λ; backend=AutoMooncake()) +#println(∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,print_times=true)) +#@btime ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol) + +#plot(plt1,plt2) + diff --git a/article/figures/time-steps.jl b/article/figures/time-steps.jl index da400c1..34e656c 100644 --- a/article/figures/time-steps.jl +++ b/article/figures/time-steps.jl @@ -63,7 +63,7 @@ n = length(funx0(λ)) sol_∂xO_flow = diagm([exp(tf*λ[1]),exp(tf*λ[2]),exp(tf*(λ[1]-λ[2]))]) sol_∂λ_flow = [λ[2]*exp(λ[1]*tf) exp(λ[1]*tf) 0. tf*exp(λ[2]*tf) - tf*exp((λ[1]-λ[2])*tf) tf*exp((λ[1]-λ[2])*tf)] + tf*exp((λ[1]-λ[2])*tf) -tf*exp((λ[1]-λ[2])*tf)] algo = Tsit5() reltol = 1.e-4; abstol = 1.e-4 @@ -96,7 +96,7 @@ sol1 = my_∂λ_flow((t0,tf),xX0,λ;reltol=RelTol,abstol=AbsTol, dt=dt) sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=reltol,abstol=abstol) dict_T_var[:var_reltol] = sol_var.t -sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=reltol,abstol=abstol,internalnorm = (u,t)->norm(u)) +sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=reltol,abstol=abstol,internalnorm = (u,t)->norm(u[:,1:1])) dict_T_var[:var_reltol_internalnorm] = sol_var.t my_Inf = prevfloat(typemax(Float64)) @@ -142,9 +142,8 @@ println("t0= ", t0) println("tf= ", tf) println("λ = ", λ) N = 5 -sol_diff_auto_flow_tf, T = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N,print_times=true) +sol_diff_auto_flow_tf, T = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, print_times=true) #, alg = Euler(),adaptive=false, dt = (tf-t0)/N,print_times=true) dict_T_auto_diff[:diff_auto_Zygote] = T -xX_auto_diff_Zygote = sol_diff_auto_flow_tf [sol_∂λ_flow sol_diff_auto_flow sol_diff_auto_flow_tf] #sort(dict_T; rev = true) diff --git a/article/figures/toto.jl b/article/figures/toto.jl new file mode 100644 index 0000000..d3f159c --- /dev/null +++ b/article/figures/toto.jl @@ -0,0 +1,21 @@ + +using DifferentiationInterface +using ForwardDiff: ForwardDiff +using Zygote: Zygote +#include("DiffFlow.jl") +#using .DiffFlow +#using OrdinaryDiffEq + +function tt(a::Vector{<:Real} ; b=1.::Real) + return a[1]+b +end + +B = [1 2 ; 3 4] +fun1(x) = B*x +t0 = 0. ; tf = 1.; x0 = [1,2]; +tspan = (t0,tf) +backend=AutoZygote() +backend1 = AutoForwardDiff() +jacobian(fun1,backend1,x0) + + diff --git a/article/main.aux b/article/main.aux index c4d354d..726a417 100644 --- a/article/main.aux +++ b/article/main.aux @@ -23,16 +23,16 @@ \@writefile{toc}{\contentsline {section}{\numberline {4}Tests on the time steps}{4}{section.4}\protected@file@percent } \@writefile{toc}{\contentsline {subsection}{\numberline {4.1}Example}{4}{subsection.4.1}\protected@file@percent } \@writefile{toc}{\contentsline {subsection}{\numberline {4.2}Step grids with variationnal equations}{4}{subsection.4.2}\protected@file@percent } -\bibstyle{plain} -\bibdata{./main} -\bibcite{Bock81}{1} \@writefile{lot}{\contentsline {table}{\numberline {1}{\ignorespaces \it Results with variationnal equations.}}{5}{table.1}\protected@file@percent } \newlabel{table:steps_var_equation}{{1}{5}{\it Results with variationnal equations}{table.1}{}} \@writefile{toc}{\contentsline {subsection}{\numberline {4.3}Step grids with automatic differentiation of the flow}{5}{subsection.4.3}\protected@file@percent } \@writefile{lot}{\contentsline {table}{\numberline {2}{\ignorespaces it Results with automatic differentiation.}}{5}{table.2}\protected@file@percent } \newlabel{table:steps_auto_diff}{{2}{5}{it Results with automatic differentiation}{table.2}{}} \@writefile{toc}{\contentsline {subsection}{\numberline {4.4}Comparaison of the solutions for the methods with the same step grid}{5}{subsection.4.4}\protected@file@percent } -\@writefile{toc}{\contentsline {subsection}{\numberline {4.5}conclusion}{5}{subsection.4.5}\protected@file@percent } +\bibstyle{plain} +\bibdata{./main} +\bibcite{Bock81}{1} \bibcite{HaNoWa93}{2} +\@writefile{toc}{\contentsline {subsection}{\numberline {4.5}conclusion}{6}{subsection.4.5}\protected@file@percent } \@writefile{toc}{\contentsline {section}{\numberline {5}Conclusion and perspectives}{6}{section.5}\protected@file@percent } \gdef \@abspage@last{6} diff --git a/article/main.log b/article/main.log index 96bcfec..a99954d 100644 --- a/article/main.log +++ b/article/main.log @@ -1,4 +1,4 @@ -This is pdfTeX, Version 3.141592653-2.6-1.40.26 (TeX Live 2024) (preloaded format=pdflatex 2024.10.11) 5 NOV 2025 23:40 +This is pdfTeX, Version 3.141592653-2.6-1.40.26 (TeX Live 2024) (preloaded format=pdflatex 2024.10.11) 6 NOV 2025 12:22 entering extended mode restricted \write18 enabled. file:line:error style messages enabled. @@ -514,7 +514,7 @@ Overfull \hbox (84.8217pt too wide) in paragraph at lines 283--283 172903, 0.30947, 0.475902, 0.666691, 0.881442, 1.0][] [] -(./main.bbl [5]) [6] (./main.aux) +[5] (./main.bbl) [6] (./main.aux) *********** LaTeX2e <2023-11-01> patch level 1 L3 programming layer <2024-02-20> @@ -533,7 +533,7 @@ Here is how much of TeX's memory you used: 32379 multiletter control sequences out of 15000+600000 568927 words of font info for 71 fonts, out of 8000000 for 9000 1141 hyphenation exceptions out of 8191 - 75i,11n,79p,541b,577s stack positions out of 10000i,1000n,20000p,200000b,200000s + 75i,14n,79p,541b,577s stack positions out of 10000i,1000n,20000p,200000b,200000s -Output written on main.pdf (6 pages, 328565 bytes). +Output written on main.pdf (6 pages, 332455 bytes). PDF statistics: 227 PDF objects out of 1000 (max. 8388607) 176 compressed objects within 2 object streams diff --git a/article/main.pdf b/article/main.pdf index b7eb3bb..a9e4727 100644 Binary files a/article/main.pdf and b/article/main.pdf differ diff --git a/article/main.synctex.gz b/article/main.synctex.gz index 7321861..cf13c27 100644 Binary files a/article/main.synctex.gz and b/article/main.synctex.gz differ diff --git a/article/main.tex b/article/main.tex index b1e263f..57f2320 100644 --- a/article/main.tex +++ b/article/main.tex @@ -320,9 +320,8 @@ \subsection{conclusion} \item The use of Zygote with the numerical integration is not good ! \end{itemize} - - - +A voir : +https://github.com/SciML/OrdinaryDiffEq.jl/blob/c174fbc1b07c252fe8ec8ad5b6e4d5fb9979c813/test/odeinterface/init_dt_vs_dopri_tests.jl#L13 %------------------------------------------------------------------------------------------------------------------ \section{Conclusion and perspectives} diff --git a/article/remarks.md b/article/remarks.md new file mode 100644 index 0000000..2702ab2 --- /dev/null +++ b/article/remarks.md @@ -0,0 +1,3 @@ +# Use of Dirrefential equation + +https://docs.sciml.ai/SciMLTutorialsOutput/html/advanced/02-advanced_ODE_solving.html \ No newline at end of file diff --git a/src/CTDiffFlow.jl b/src/CTDiffFlow.jl index 6d28083..2929288 100644 --- a/src/CTDiffFlow.jl +++ b/src/CTDiffFlow.jl @@ -182,14 +182,11 @@ function build_∂λ_flow(rhs::Function,t0::Real,x0::Function,tf::Real, λ::Vect function _flow(λ) ivp = ODEProblem(rhs, x0(λ), (t0,tf), λ) algo = get(ode_kwargs, :alg, Tsit5()) - println("algo = ", algo) sol = solve(ivp, alg=algo; ode_kwargs...) T = sol.t - println(sol.u[end]) return sol.u[end] end if print_times - println("backend = ", backend) return jacobian(_flow,backend,λ0), T else return jacobian(_flow,backend,λ0) diff --git a/src/myode43/figHairer.jl b/src/myode43/figHairer.jl new file mode 100644 index 0000000..d699483 --- /dev/null +++ b/src/myode43/figHairer.jl @@ -0,0 +1,94 @@ +# +# ~gergaud/ENS/Control/ODE/pasvar/figHairer.m +# +# Auteurs: Joseph GERGAUD +# Date: fev. 2007 +# Adresse: INP-ENSEEIHT-IRIT-UMR CNRS 5055 +# 2, rue Camichel 31071 Toulouse FRANCE +# Email: gergaud@enseeiht.fr +#*************************************************************************** +# +# pour visualiser les graphiques du livre +# ref: Hairer page 170 T1 +# Nombre de pas acceptes et rejetes pas variables +include("myode43_for_graph.jl") +using Plots + +function bruss(x,par,t) + λ = par[1] + x₁ = x[1] + x₂ = x[2] + return [1+x₁^2*x₂-(λ+1)*x₁ , λ*x₁-x₁^2*x₂] +end + +t0=0;tf=20; +# +Atol=1.e-4; Rtol=1.e-4; +x0=[1.5,3]; +par = [3.] +T,X,haccept,Trej,hrej,Ttot,err_loc_est,err_loc_exact,err_abs = myode43(bruss,x0,par,[t0, tf],Rtol,Atol) + +N = length(T) +println("N = ", N) +println(size(X)) + +X1 = [X[i][1] for i in 1:N] +p1 = plot(T,X1,labels=false) + +scatter!(p1,T,X1,c=:blue,labels=false) +X2 = [X[i][2] for i in 1:N] +plot!(p1,T,X2,c=:green,labels=false) +scatter!(p1,T,X2,c=:green,labels=false,ylabel="x_1, x_2") +# + +npasaccept = length(haccept); +p2 = plot(T[1:npasaccept],haccept,yscale = :log10,labels=false); +scatter!(p2,T[1:npasaccept],haccept,yscale = :log10,c=:blue,labels="accepted steps"); +scatter!(p2,Trej,hrej,yscale = :log10,c=:red,labels="reject steps"); + + + +p3 = plot(Ttot[2:end],err_loc_est,yscale = :log10,c=:blue,labels="local estimated errors") +scatter!(p3,Ttot[2:end],err_loc_est,yscale = :log10,c=:blue,labels=false, +markersize = 3, markershape = :x) +plot!(p3,[t0, tf],[Atol, Atol],c=:black,labels=false) + +plot!(p3,Ttot[2:end],err_loc_exact,c=:green,yscale = :log10,labels="local estimated errors") +scatter!(p3,Ttot[2:end],err_loc_exact,c=:green,yscale = :log10,labels=false, +markersize = 3,markershape = :x) +plot!(p3,Ttot[2:end],err_abs,c=:black,yscale = :log10,labels="global exact errors") +scatter!(p3,Ttot[2:end],err_abs,c=:black,yscale = :log10,labels=false,markershape = :x, +markersize = 3) + +plot(p1,p2,p3,layout=(3,1)) + +#= +# +# Equation avec discintinuit� du second membre +# Figure page 197 avec mon programme pas variable +t0=0;tf=1; +# +Atol=0.7e-3; Rtol=0.7e-3; +y0=0.3; +[T,Y,haccept,Trej,hrej,Ttot,err_loc_est,err_loc_exact,err_abs]=myode43(@phi_disc,[t0 tf],y0,Rtol,Atol); +figure; +j = 1; +for i=1:length(Ttot), + if T(j)==Ttot(i), + I(i)=1; + j=j+1; + else + I(i)=0; + end; +end; +h1=plot(Ttot(find(I==1)),find(I==1)-1,'o') +hold on +h2=plot(Ttot(find(I==0)),find(I==0)-1,'xr') +plot(Ttot,0:length(Ttot)-1) +xlabel('t'); +ylabel('Nombre de pas'); +legend([h1 h2],'pas accept�s','pas rejet�s') +print('fig_disc1','-depsc') + + +=# \ No newline at end of file diff --git a/src/myode43/myode43.jl b/src/myode43/myode43.jl new file mode 100644 index 0000000..1ad61db --- /dev/null +++ b/src/myode43/myode43.jl @@ -0,0 +1,112 @@ +# +# Numerical integration with control step size +# algorithm : see Hairer Tome 1 page +# +# +using LinearAlgebra +function inith(rhs,x0,par,t0,Atol,Rtol) + # + # + n = length(x0); + sc = Atol .+ abs.(x0) .* Rtol; + println("sc = ", sc) + println("x0 = ", x0) + k0 = rhs(x0, par, t0); + println(x0./sc) + println(norm(x0./sc)) + d0 = norm(x0./sc)/sqrt(n); # normalement c'est la bonne valeur cf page 168 + d1 = norm(k0./sc)/sqrt(n); # normalement c'est la bonne valeur cf page 168 + #d0 = norm(x0./sc); + #d1 = norm(k0./sc); + h0 = 0.01*(d0/d1); + if (d0 < 1.e-5) || (d1 < 1.e-5) + h0 = 1.e-6; + end; + x1 = x0 + h0*k0; + k1 = rhs(x1,par,t0+h0); + d2 = norm((k1-k0)./sc)/sqrt(n)/h0; # normalement c'est la bonne valeur cf page 168 + #d2 = norm((k1-k0)./sc)/h0; + if max(d1,d2) < 1.e-15 + h1 = max(1.e-6,h0*1.e-3); + else + h1 = (0.01/max(d1,d2))^(1/4); + end + h = min(100*h0,h1); + return h +end + + +function myode43(rhs,x0,par,t0tf,Rtol,Atol) +# +# Initialisation + Scal_Type=eltype(par) + p=4; # ordre + t0=t0tf[1]; tf=t0tf[2]; + hmax=tf-t0; + Npasmax=1000; + n = length(x0); + T = [Scal_Type(t0)]; + #T = [t0]; + X = [Scal_Type.(x0)] + # + # + # step initialisation + #h=0.03; + h=inith(rhs,x0,par,t0,Atol,Rtol) + #println("h= ",h) + # + nstep=0; + t=t0 + x=x0; fin=0; + while (nstep < Npasmax) && (fin == 0) + nstep = nstep+1 + # Runge-Kutta method + k1 = rhs(x,par,t); + k2 = rhs(x+(h/3)*k1,par,t+h/3,); + k3 = rhs(x+h*(-k1/3+k2),par,t+2*h/3); + k4 = rhs(x+h*(k1-k2+k3),par,t+h); + x1 = x+(h/8)*(k1+3*k2+3*k3+k4); + xhat1 = x+(h/12)*(k1+6*k2+3*k3+2*rhs(x1,par,t)); + # err + sc = Atol .+ max.(abs.(x),abs.(x1)) .* Rtol; + #println("sc = ", sc) + err = norm((x1-xhat1)./sc)/sqrt(n); + # println("err = ", err) + # + # calcul du pas + if (err < 1) + t = t+h; x = x1; + #println("t= ",t) + #println("T= ",T) + #println("X= ",X) + push!(T, t) + #println("T= ",T) + #println("x1= ",x1) + push!(X, x1) + #println("X= ",X) + if (t > tf-h/2) + fin = 1; + end; + end; + h = h*min(5,max(0.2,0.9*(1/err)^(1/p))); + if (t+h > tf) + h = tf-t; + end; + end; + + return T,X + end + + function get_Xij(X,i,j) + """ + Get from the vector of vector X the i,j composante for all times + """ + N = length(X) + Xij = zeros(N) + for l in 1:N + Xij[l] = X[l][i,j] + end + return Xij + end + + diff --git a/src/myode43/myode43_for_graph.jl b/src/myode43/myode43_for_graph.jl new file mode 100644 index 0000000..11cf063 --- /dev/null +++ b/src/myode43/myode43_for_graph.jl @@ -0,0 +1,109 @@ +# +# Numerical integration with control step size +# algorithm : see Hairer Tome 1 page +# +# +using LinearAlgebra +function inith(rhs,x0,par,t0,Atol,Rtol) + # + # + n = length(x0); + sc = Atol .+ abs.(x0) .* Rtol; + k0 = rhs(x0, par, t0); + #d0 = norm(x0./sc)/sqrt(n); # normalement c'est la bonne valeur cf page 168 + #d1 = norm(k0./sc)/sqrt(n); # normalement c'est la bonne valeur cf page 168 + d0 = norm(x0./sc); + d1 = norm(k0./sc); + h0 = 0.01*(d0/d1); + if (d0 < 1.e-5) || (d1 < 1.e-5) + h0 = 1.e-6; + end; + x1 = x0 + h0*k0; + k1 = rhs(x1,par,t0+h0); + # d2 = norm((k1-k0)./sc)/sqrt(n)/h0; # normalement c'est la bonne valeur cf page 168 + d2 = norm((k1-k0)./sc)/h0; + if max(d1,d2) < 1.e-15 + h1 = max(1.e-6,h0*1.e-3); + else + h1 = (0.01/max(d1,d2))^(1/4); + end + h = min(100*h0,h1); + return h +end + + +function myode43(rhs,x0,par,t0tf,Rtol,Atol) +# +# Initialisation + p=4; # ordre + t0=t0tf[1]; tf=t0tf[2]; + hmax=tf-t0; + Npasmax=1000; + n=length(x0); + naccept=0; nrej=0; + T=Float64[t0]; X=[x0]; haccept=Float64[]; + Trej=Float64[]; hrej=Float64[]; + Ttot=Float64[t0]; err_loc_est=Float64[]; err_loc_exact=Float64[]; err_abs=Float64[]; + # + # + # step initialisation + #h=0.03; + h=inith(rhs,x0,par,t0,Atol,Rtol) + # initialisation pour avoir les erreurs ''exactes'' + RelTol = 1.e-7 + AbsTol = 1.e-7 + # + nstep=0; t=t0 + x=x0; fin=0; + while (nstep < Npasmax) && (fin == 0) + nstep = nstep+1 + # Runge-Kutta method + k1 = rhs(x,par,t); + k2 = rhs(x+(h/3)*k1,par,t+h/3,); + k3 = rhs(x+h*(-k1/3+k2),par,t+2*h/3); + k4 = rhs(x+h*(k1-k2+k3),par,t+h); + x1 = x+(h/8)*(k1+3*k2+3*k3+k4); + xhat1 = x+(h/12)*(k1+6*k2+3*k3+2*rhs(x1,par,t)); + # err + sc = Atol .+ max.(abs.(x),abs.(x1)) .* Rtol; + #println("sc = ", sc) + err = norm((x1-xhat1)./sc)/sqrt(n); + #println("err = ", err) + # err_loc_est = [err_loc_est ; norm(y1-yhat1)]; + push!(err_loc_est, err*Atol) + push!(Ttot, t+h) + algo = Tsit5() + ivp = ODEProblem(rhs, x, (t,t+h), par) + sol = solve(ivp, alg=algo, reltol = RelTol, abstol = AbsTol) + err_exact = norm((x1-sol.u[end])./sc)/sqrt(n); + push!(err_loc_exact, err_exact*Atol) + ivp = ODEProblem(rhs, x0, (t0,t+h), par) + sol = solve(ivp, alg=algo, reltol = RelTol, abstol = AbsTol) + #[TT,YY] = ode45(f,[t0 t+h],y0,optionode); + err_globale = norm((x1-sol.u[end])./sc)/sqrt(n); + push!(err_abs, err_globale*Atol) + # + # calcul du pas + if (err < 1) + naccept = 1+naccept; + t = t+h; x = x1; + push!(haccept, h) + push!(T, t) + push!(X, x1) + if (t > tf-h/2) + fin = 1; + end; + else + nrej = nrej+1; + push!(Trej,t+h) + push!(hrej,h) + end; + h = h*min(5,max(0.2,0.9*(1/err)^(1/p))); + if (t+h > tf) + h = tf-t; + end; + end; + + return T,X,haccept,Trej,hrej,Ttot,err_loc_est,err_loc_exact,err_abs + end +