diff --git a/article/figures/chris.jl b/article/figures/chris.jl index 1b062db..d4a4fc1 100644 --- a/article/figures/chris.jl +++ b/article/figures/chris.jl @@ -32,25 +32,30 @@ function hinit(F, x0, t0, tend, p, reltol, abstol) 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 -A(λ) = [λ[1] 0 0 ; 0 λ[2] 0 ; 0 0 λ[1]-λ[2]] -fun1(t,x) = A(λ)*x +function 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]]] +end -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.; +t0 = 0.0; +tend = 2.0; p = 2 -x0 = [0.,1.,1] -reltol = 1.e-4; abstol = 1.e-4 +x0 = [0.0, 1.0, 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]] +xX0 = [funx0(λ) [0 1; 0 0; 0 0]] hinit(fun2, xX0, t0, tend, p, reltol, abstol) -RelTol = reltol*ones(n,p+1) +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) +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 index 181364c..fea9733 100644 --- a/article/figures/test-dopri-DP5-julia.jl +++ b/article/figures/test-dopri-DP5-julia.jl @@ -1,5 +1,4 @@ -using OrdinaryDiffEq, Test, DiffEqDevTools, - ODEInterface, ODEInterfaceDiffEq +using OrdinaryDiffEq, Test, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq import ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear @@ -20,7 +19,7 @@ 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))) +sol = solve(prob, DP5(); internalnorm=(u, t) -> sqrt(sum(abs2, u))) # Change the norm due to error in dopri5.f sol2 = solve(prob, dopri5()) @@ -40,9 +39,9 @@ println("sol.t[2] - sol2.t[2] = ", 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))) +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 +@test sol.t[2] ≈ sol2.t[2] diff --git a/article/figures/test_myode43.jl b/article/figures/test_myode43.jl index 5dece0a..cb2b624 100644 --- a/article/figures/test_myode43.jl +++ b/article/figures/test_myode43.jl @@ -1,6 +1,5 @@ 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₀(λ) @@ -9,48 +8,52 @@ include("../../src/myode43/myode43.jl") # # ∂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]; +A(λ) = [λ[1] 0 0; 0 λ[2] 0; 0 0 λ[1]-λ[2]] +fun1(x, λ, t) = A(λ)*x + +function 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]]] +end + +t0 = 0.0; +tf = 2.0; +tspan = (t0, tf); +λ = [1.0, 2]; +p = length(λ); +funx0(λ) = [λ[2], 1.0, 1]; x0 = funx0(λ) n = length(x0) -reltol = 1.e-4; abstol = 1.e-4 -myode43(fun1,x0,λ,tspan,reltol,abstol) +reltol = 1.e-4; +abstol = 1.e-4 +myode43(fun1, x0, λ, tspan, reltol, abstol) -xX0 = [x0 [0 1 ; 0 0 ; 0 0]] +xX0 = [x0 [0 1; 0 0; 0 0]] -myode43(fun2,xX0,λ,tspan,reltol,abstol) +myode43(fun2, xX0, λ, tspan, reltol, abstol) -RelTol = reltol*ones(n,p+1) +RelTol = reltol*ones(n, p+1) AbsTol = RelTol -myode43(fun2,xX0,λ,tspan,RelTol,AbsTol) +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) +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(fun1, x0, λ, t0, abstol, reltol) -inith(fun2,xX0,λ,t0,abstol,reltol) -RelTol = reltol*ones(n,p+1) +inith(fun2, xX0, λ, t0, abstol, reltol) +RelTol = reltol*ones(n, p+1) AbsTol = RelTol -inith(fun2,xX0,λ,t0,AbsTol,RelTol) - +inith(fun2, xX0, λ, t0, AbsTol, RelTol) epsi = eps() epsi = 0 -xX0 = [x0 [epsi 1 ; epsi epsi ; epsi epsi]] +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) - +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 index d9080a3..875c398 100644 --- a/article/figures/test_zygote.jl +++ b/article/figures/test_zygote.jl @@ -13,81 +13,120 @@ include("../../src/CTDiffFlow.jl") using .CTDiffFlow include("../../src/myode43/myode43.jl") -function my_euler(fun,t0,x0,tf,λ,N) - ti = t0; xi = x0 +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) + 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]; +A(λ) = [λ[1] 0 0; 0 λ[2] 0; 0 0 λ[1]-λ[2]] +fun1(x, λ, t) = A(λ)*x +funx0(λ) = [λ[2], 1.0, 1]; N = 10 -t0 = 0. ; tf = 1.; tspan = (t0,tf); -λ = [1.,2]; +t0 = 0.0; +tf = 1.0; +tspan = (t0, tf); +λ = [1.0, 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)] - +sol_∂λ_flow = [ + λ[2]*exp(λ[1]*tf) exp(λ[1]*tf) + 0.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}[]) +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) +_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(),λ))) +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] + 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(),λ))) +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()) +∂λ_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))) +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()) +∂λ_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))) +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) + 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(),λ))) +push!(df, ("ForwardDiff", "my-ode43", jacobian(_flow, AutoForwardDiff(), λ))) -jacobian(_flow,AutoForwardDiff(),λ) +jacobian(_flow, AutoForwardDiff(), λ) -latexify(df,arraystyle=:pmatrix,env=:tabular) +latexify(df; arraystyle=:pmatrix, env=:tabular) -sol_Zygote = df[in.(df.AutoDiff, Ref(["solution","Zygote"])),:] -sol_ForwardDiff = df[in.(df.AutoDiff, Ref(["solution", "ForwardDiff"])),:] +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 index f9b7684..02cefc5 100644 --- a/article/figures/time-steps-bruss.jl +++ b/article/figures/time-steps-bruss.jl @@ -15,29 +15,32 @@ 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₂] +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.] +t0 = 0.0; +tf = 1.0 +tspan = (t0, tf) +λ = [3.0] p = length(λ) x01 = 1.3 funx0(λ) = [x01, λ[1]] tol = 1.e-4 n = 2 -plt1 = plot(); plt2 = plot() +plt1 = plot(); +plt2 = plot() algo = Tsit5() -reltol = 1.e-4; abstol = 1.e-4 +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) +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) @@ -48,18 +51,20 @@ dt = sol_ivp.t[2] println("Vatiational equation") # ------------------------------ -∂λ_flow_var = CTDiffFlow.build_∂λ_flow_var(bruss,t0,funx0,tf,λ) +∂λ_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)) +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) ) +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) @@ -70,8 +75,10 @@ dict_T[:var_RelTol] = sol_var.t 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) +∂λ_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) @@ -82,8 +89,10 @@ println(∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,internalnorm = (u, =# 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) +∂λ_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) @@ -99,4 +108,3 @@ println("With Mooncake") #@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 34e656c..29363b4 100644 --- a/article/figures/time-steps.jl +++ b/article/figures/time-steps.jl @@ -24,54 +24,60 @@ using .CTDiffFlow # # ∂x0_flow -A(λ) = [λ[1] 0 0 ; 0 λ[2] 0 ; 0 0 λ[1]-λ[2]] -fun1(x,λ,t) = A(λ)*x -function my_∂x0_flow(tspan,x0,λ; ode_kwargs...) +A(λ) = [λ[1] 0 0; 0 λ[2] 0; 0 0 λ[1]-λ[2]] +fun1(x, λ, t) = A(λ)*x +function my_∂x0_flow(tspan, x0, λ; ode_kwargs...) n = length(x0) In = Matrix(I(n)) - ivp_fun1 = ODEProblem(fun1, In , tspan, λ) + ivp_fun1 = ODEProblem(fun1, In, tspan, λ) alg = get(ode_kwargs, :alg, Tsit5()) - sol = solve(ivp_fun1, alg=alg; ode_kwargs...) + sol = solve(ivp_fun1; alg=alg, ode_kwargs...) return sol end -function my_∂x0_flow(t0,x0,tf,λ; ode_kwargs...) - sol = my_∂x0_flow((t0,tf),x0,λ; ode_kwargs...) +function my_∂x0_flow(t0, x0, tf, λ; ode_kwargs...) + sol = my_∂x0_flow((t0, tf), x0, λ; ode_kwargs...) return sol.u[end] end -function my_∂λ_flow(tspan,xX0,λ; ode_kwargs...) +function my_∂λ_flow(tspan, xX0, λ; ode_kwargs...) #n,p = size(xX0) - 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]] ] - - ivp_fun = ODEProblem(fun1, xX0 , tspan, λ) + 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]]] + + ivp_fun = ODEProblem(fun1, xX0, tspan, λ) alg = get(ode_kwargs, :alg, Tsit5()) - sol = solve(ivp_fun, alg=alg; ode_kwargs...) + sol = solve(ivp_fun; alg=alg, ode_kwargs...) return sol end -function my_∂λ_flow(t0,xX0,tf,λ; ode_kwargs...) - sol = my_∂λ_flow((t0,tf),xX0,λ; ode_kwargs...) +function my_∂λ_flow(t0, xX0, tf, λ; ode_kwargs...) + sol = my_∂λ_flow((t0, tf), xX0, λ; ode_kwargs...) return sol.u[end] end - -λ = [1.,2]; p = length(λ); -t0 = 0. ; tf = 1.; tspan = (t0,tf); -funx0(λ) = [λ[2],1.,1]; +λ = [1.0, 2]; +p = length(λ); +t0 = 0.0; +tf = 1.0; +tspan = (t0, tf); +funx0(λ) = [λ[2], 1.0, 1]; 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)] +sol_∂xO_flow = diagm([exp(tf*λ[1]), exp(tf*λ[2]), exp(tf*(λ[1]-λ[2]))]) +sol_∂λ_flow = [ + λ[2]*exp(λ[1]*tf) exp(λ[1]*tf) + 0.0 tf*exp(λ[2]*tf) + tf*exp((λ[1]-λ[2])*tf) -tf*exp((λ[1]-λ[2])*tf) +] algo = Tsit5() reltol = 1.e-4; abstol = 1.e-4 # Flow #x0λ = funx0(λ) -ivp = ODEProblem(fun1, funx0(λ), (t0,tf), λ) -sol_ivp = solve(ivp, alg=algo; reltol=reltol,abstol=abstol) +ivp = ODEProblem(fun1, 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) @@ -83,66 +89,79 @@ dt = sol_ivp.t[2] println("Vatiational equation") # ------------------------------ -xX0 = [funx0(λ) [0 1 ; 0 0 ; 0 0]] -println(my_∂λ_flow(t0,xX0,tf,λ)) +xX0 = [funx0(λ) [0 1; 0 0; 0 0]] +println(my_∂λ_flow(t0, xX0, tf, λ)) -RelTol = [reltol*ones(n,1) Inf*ones(n,2)]/sqrt(p+1) -AbsTol = [abstol*ones(n,1) Inf*ones(n,2)]/sqrt(p+1) +RelTol = [reltol*ones(n, 1) Inf*ones(n, 2)]/sqrt(p+1) +AbsTol = [abstol*ones(n, 1) Inf*ones(n, 2)]/sqrt(p+1) # plante car 0.*Inf = NaN lors de l'initialisation du pas -sol1 = my_∂λ_flow((t0,tf),xX0,λ;reltol=RelTol,abstol=AbsTol, dt=dt) +sol1 = my_∂λ_flow((t0, tf), xX0, λ; reltol=RelTol, abstol=AbsTol, dt=dt) -∂λ_flow_var = CTDiffFlow.build_∂λ_flow_var(fun1,t0,funx0,tf,λ) +∂λ_flow_var = CTDiffFlow.build_∂λ_flow_var(fun1, t0, funx0, tf, λ) -sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=reltol,abstol=abstol) +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[:,1:1])) +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)) -RelTol = [reltol*ones(n,1) my_Inf*ones(n,2)]/sqrt(p+1) -AbsTol = [abstol*ones(n,1) my_Inf*ones(n,2)]/sqrt(p+1) - +RelTol = [reltol*ones(n, 1) my_Inf*ones(n, 2)]/sqrt(p+1) +AbsTol = [abstol*ones(n, 1) my_Inf*ones(n, 2)]/sqrt(p+1) println("AbsTol = ", AbsTol) println("RelTol = ", RelTol) -sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=RelTol,abstol=AbsTol) +sol_var = ∂λ_flow_var((t0, tf), funx0, λ; reltol=RelTol, abstol=AbsTol) println(sol_ivp.t - sol_var.t) #println("dt = ",dt) -sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=RelTol,abstol=AbsTol,dt=dt) +sol_var = ∂λ_flow_var((t0, tf), funx0, λ; reltol=RelTol, abstol=AbsTol, dt=dt) println("RelTol = ", RelTol) println("Times = ", sol_var.t) println(sol_ivp.t - sol_var.t) dict_T_var[:var_RelTol] = sol_var.t -xX_var_tf = sol_var.u[end][:,2:end] - +xX_var_tf = sol_var.u[end][:, 2:end] # ne fonctionne pas #@btime ∂λ_flow_var(t0,funx0,tf,λ;reltol=reltol,abstol=abstol) println("Automatic differentiation") println("with ForwardDiff") -∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1,t0,funx0,tf,λ) -sol_diff_auto_flow, T = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,print_times=true) +∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1, t0, funx0, tf, λ) +sol_diff_auto_flow, T = ∂λ_flow( + t0, funx0, tf, λ; reltol=reltol, abstol=abstol, print_times=true +) println("Times : ", T) dict_T_auto_diff[: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)") -sol_diff_auto_flow, T = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,internalnorm = (u,t)->norm(u),print_times=true) +sol_diff_auto_flow, T = ∂λ_flow( + t0, + funx0, + tf, + λ; + reltol=reltol, + abstol=abstol, + internalnorm=(u, t)->norm(u), + print_times=true, +) println("Times : ", T) dict_T_auto_diff[:diff_auto_ForwardDiff_internalnorm] = T #@btime ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol) println("With Zygote") # with Zygote -∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1,t0,funx0,tf,λ; backend=AutoZygote()) +∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1, t0, funx0, tf, λ; backend=AutoZygote()) println("t0= ", t0) println("tf= ", tf) println("λ = ", λ) N = 5 -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) +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 [sol_∂λ_flow sol_diff_auto_flow sol_diff_auto_flow_tf] @@ -159,6 +178,3 @@ println("With Mooncake") #∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1,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) - - - diff --git a/article/figures/toto.jl b/article/figures/toto.jl index d3f159c..6a9bf47 100644 --- a/article/figures/toto.jl +++ b/article/figures/toto.jl @@ -6,16 +6,16 @@ using Zygote: Zygote #using .DiffFlow #using OrdinaryDiffEq -function tt(a::Vector{<:Real} ; b=1.::Real) - return a[1]+b +function tt(a::Vector{<:Real}; b=1.0::Real) + return a[1]+b end -B = [1 2 ; 3 4] +B = [1 2; 3 4] fun1(x) = B*x -t0 = 0. ; tf = 1.; x0 = [1,2]; -tspan = (t0,tf) +t0 = 0.0; +tf = 1.0; +x0 = [1, 2]; +tspan = (t0, tf) backend=AutoZygote() backend1 = AutoForwardDiff() -jacobian(fun1,backend1,x0) - - +jacobian(fun1, backend1, x0) diff --git a/formatter.lock b/formatter.lock new file mode 100644 index 0000000..77f0d0c --- /dev/null +++ b/formatter.lock @@ -0,0 +1 @@ +JuliaFormatter-1 diff --git a/src/CTDiffFlow.jl b/src/CTDiffFlow.jl index 2929288..ab887e2 100644 --- a/src/CTDiffFlow.jl +++ b/src/CTDiffFlow.jl @@ -4,7 +4,7 @@ using DifferentiationInterface using OrdinaryDiffEq using LinearAlgebra -function built_rhs_var(rhs::Function; wrt = :x0, backend = AutoForwardDiff() ) +function built_rhs_var(rhs::Function; wrt=:x0, backend=AutoForwardDiff()) """ Built the second member of the variational equations input @@ -25,61 +25,83 @@ function built_rhs_var(rhs::Function; wrt = :x0, backend = AutoForwardDiff() ) The 2:end column is the state X of the variational equations """ - @assert (wrt == :x0 || wrt == :λ || wrt ==:t0) "Error the wrt optional argument of the built_rhs_var function is not equal to :x0 or :λ ot :t0" - fun_x(x,λ,t) = jacobian(x -> rhs(x,λ,t), backend, x) + @assert (wrt == :x0 || wrt == :λ || wrt == :t0) "Error the wrt optional argument of the built_rhs_var function is not equal to :x0 or :λ ot :t0" + fun_x(x, λ, t) = jacobian(x -> rhs(x, λ, t), backend, x) if wrt == :λ - fun_λ(x,λ,t) = jacobian(λ -> rhs(x,λ,t), backend, λ) - function rhs_var_λ(xX,λ,t) - x = xX[:,1] - X = xX[:,2:end] - xpoint = rhs(x,λ,t) - Xpoint = fun_x(x,λ,t)*X + fun_λ(x,λ,t) + fun_λ(x, λ, t) = jacobian(λ -> rhs(x, λ, t), backend, λ) + function rhs_var_λ(xX, λ, t) + x = xX[:, 1] + X = xX[:, 2:end] + xpoint = rhs(x, λ, t) + Xpoint = fun_x(x, λ, t)*X + fun_λ(x, λ, t) return [xpoint Xpoint] end return rhs_var_λ else - function rhs_var_x0(xX,λ,t) - x = xX[:,1] - X = xX[:,2:end] - xpoint = rhs(x,λ,t) - Xpoint = fun_x(x,λ,t)*X + function rhs_var_x0(xX, λ, t) + x = xX[:, 1] + X = xX[:, 2:end] + xpoint = rhs(x, λ, t) + Xpoint = fun_x(x, λ, t)*X return [xpoint Xpoint] end - return rhs_var_x0 + return rhs_var_x0 end end - # Integration of the variational equations # derivatives with respect to x0 -function build_∂x0_flow_var(rhs::Function,t0::Real,x0::Vector{<:Real},tf::Real, λ::Vector{<:Real}; backend = AutoForwardDiff()) #,print_step=false) - - rhs_var = built_rhs_var(rhs , wrt = :x0, backend = backend) +function build_∂x0_flow_var( + rhs::Function, + t0::Real, + x0::Vector{<:Real}, + tf::Real, + λ::Vector{<:Real}; + backend=AutoForwardDiff(), +) #,print_step=false) + rhs_var = built_rhs_var(rhs; wrt=:x0, backend=backend) - function ∂x0_flow(tspan::Tuple{<:Real,<:Real},x0::Vector{<:Real}, λ::Vector{<:Real}; print_times=false, ode_kwargs...) + function ∂x0_flow( + tspan::Tuple{<:Real,<:Real}, + x0::Vector{<:Real}, + λ::Vector{<:Real}; + print_times=false, + ode_kwargs..., + ) n = length(x0) x0δx0 = [x0 Matrix(I(n))] ivp = ODEProblem(rhs_var, x0δx0, tspan, λ) algo = get(ode_kwargs, :alg, Tsit5()) - sol = solve(ivp, alg=algo; ode_kwargs...) + sol = solve(ivp; alg=algo, ode_kwargs...) return sol end - function ∂x0_flow(t0::Real,x0::Vector{<:Real}, tf::Real, λ::Vector{<:Real}; ode_kwargs...) + function ∂x0_flow( + t0::Real, x0::Vector{<:Real}, tf::Real, λ::Vector{<:Real}; ode_kwargs... + ) #println("reltol = ", get(ode_kwargs, :reltol, 1.e-3)) #println("abstol = ", get(ode_kwargs, :abstol, 1.e-6)) - sol = ∂x0_flow((t0,tf),x0,λ; ode_kwargs...) - return sol.u[end][:,2:end] - end + sol = ∂x0_flow((t0, tf), x0, λ; ode_kwargs...) + return sol.u[end][:, 2:end] + end return ∂x0_flow end # derivatives with respect to λ -function build_∂λ_flow_var(rhs::Function,t0::Real,x0::Function,tf::Real, λ::Vector{<:Real}; backend = AutoForwardDiff()) #,print_step=false) - Jacλx0 = jacobian(x0,backend,λ) - rhs_var = built_rhs_var(rhs , wrt = :λ, backend = backend) +function build_∂λ_flow_var( + rhs::Function, + t0::Real, + x0::Function, + tf::Real, + λ::Vector{<:Real}; + backend=AutoForwardDiff(), +) #,print_step=false) + Jacλx0 = jacobian(x0, backend, λ) + rhs_var = built_rhs_var(rhs; wrt=:λ, backend=backend) - function ∂λ_flow(tspan::Tuple{<:Real,<:Real},x0::Function, λ::Vector{<:Real}; ode_kwargs...) + function ∂λ_flow( + tspan::Tuple{<:Real,<:Real}, x0::Function, λ::Vector{<:Real}; ode_kwargs... + ) x0λ = x0(λ) #= n = length(x0λ); p = length(λ); @@ -87,22 +109,29 @@ function build_∂λ_flow_var(rhs::Function,t0::Real,x0::Function,tf::Real, λ:: x0δλ0 = [x0λ Jacλx0] ivp = ODEProblem(rhs_var, x0δλ0, tspan, λ) algo = get(ode_kwargs, :alg, Tsit5()) - sol = solve(ivp, alg=algo; ode_kwargs...) + sol = solve(ivp; alg=algo, ode_kwargs...) return sol end - function ∂λ_flow(t0::Real,x0::Function, tf::Real, λ::Vector{<:Real}; ode_kwargs...) + function ∂λ_flow(t0::Real, x0::Function, tf::Real, λ::Vector{<:Real}; ode_kwargs...) #println("reltol = ", get(ode_kwargs, :reltol, 1.e-3)) #println("abstol = ", get(ode_kwargs, :abstol, 1.e-6)) - sol = ∂λ_flow((t0,tf),x0,λ; ode_kwargs...) - return sol.u[end][:,2:end] - end + sol = ∂λ_flow((t0, tf), x0, λ; ode_kwargs...) + return sol.u[end][:, 2:end] + end return ∂λ_flow end # Automatic differentiation on the flow # derivatives with respect to x0 -function build_∂x0_flow(rhs::Function,t0::Real,x0::Vector{<:Real},tf::Real, λ::Vector{<:Real}; backend = AutoForwardDiff()) +function build_∂x0_flow( + rhs::Function, + t0::Real, + x0::Vector{<:Real}, + tf::Real, + λ::Vector{<:Real}; + backend=AutoForwardDiff(), +) #= function ∂x0_flow(tspan::Tuple{<:Real,<:Real},x0::Vector{<:Real}, δx0::Vector{<:Real}, λ::Vector{<:Real}; ode_kwargs...) @@ -115,84 +144,116 @@ function build_∂x0_flow(rhs::Function,t0::Real,x0::Vector{<:Real},tf::Real, λ return jacobian(flow,backend,x0) end =# -#= - function ∂x0_flow(tspan::Tuple{<:Real,<:Real},x0::Vector{<:Real}, δx0::Matrix{<:Real}, λ::Vector{<:Real}; ode_kwargs...) - function flow(x0) - ivp = ODEProblem(rhs, x0, tspan, λ) - sol = solve(ivp, alg=algo; ode_kwargs...) - return sol.u[end] + #= + function ∂x0_flow(tspan::Tuple{<:Real,<:Real},x0::Vector{<:Real}, δx0::Matrix{<:Real}, λ::Vector{<:Real}; ode_kwargs...) + function flow(x0) + ivp = ODEProblem(rhs, x0, tspan, λ) + sol = solve(ivp, alg=algo; ode_kwargs...) + return sol.u[end] + end + return jacobian(flow,backend,x0)*δx0 end - return jacobian(flow,backend,x0)*δx0 - end - =# + =# - - function ∂x0_flow(t0::Real,x0::Vector{<:Real}, tf::Real, λ::Vector{<:Real}; print_times=false, ode_kwargs...) + function ∂x0_flow( + t0::Real, + x0::Vector{<:Real}, + tf::Real, + λ::Vector{<:Real}; + print_times=false, + ode_kwargs..., + ) function _flow(x0) - ivp = ODEProblem(rhs, x0, (t0,tf), λ) + ivp = ODEProblem(rhs, x0, (t0, tf), λ) algo = get(ode_kwargs, :alg, Tsit5()) - sol = solve(ivp, alg=algo; ode_kwargs...) + sol = solve(ivp; alg=algo, ode_kwargs...) return sol.u[end] end if print_times - return jacobian(_flow,backend,x0), sol.t + return jacobian(_flow, backend, x0), sol.t else - return jacobian(_flow,backend,x0) + return jacobian(_flow, backend, x0) end - end - - function ∂x0_flow(t0::Real,x0::Vector{<:Real}, tf::Real, λ::Vector{<:Real}; print_times=false, ode_kwargs...) + end + + function ∂x0_flow( + t0::Real, + x0::Vector{<:Real}, + tf::Real, + λ::Vector{<:Real}; + print_times=false, + ode_kwargs..., + ) function _flow(x0) - ivp = ODEProblem(rhs, x0, (t0,tf), λ) + ivp = ODEProblem(rhs, x0, (t0, tf), λ) algo = get(ode_kwargs, :alg, Tsit5()) - sol = solve(ivp, alg=algo; ode_kwargs...) + sol = solve(ivp; alg=algo, ode_kwargs...) return sol.u[end] end if print_times - return jacobian(_flow,backend,x0), sol.t + return jacobian(_flow, backend, x0), sol.t else - return jacobian(_flow,backend,x0) + return jacobian(_flow, backend, x0) end - end + end ∂x0_flow end # derivatives with respect to λ -function build_∂λ_flow(rhs::Function,t0::Real,x0::Function,tf::Real, λ::Vector{<:Real}; backend = AutoForwardDiff()) - - function ∂λ_flow(tspan::Tuple{<:Real,<:Real},x0::Function, λ::Vector{<:Real}; print_times=false, ode_kwargs...) +function build_∂λ_flow( + rhs::Function, + t0::Real, + x0::Function, + tf::Real, + λ::Vector{<:Real}; + backend=AutoForwardDiff(), +) + function ∂λ_flow( + tspan::Tuple{<:Real,<:Real}, + x0::Function, + λ::Vector{<:Real}; + print_times=false, + ode_kwargs..., + ) T = [] function _flow(λ) ivp = ODEProblem(rhs, x0(λ), tspan, λ) algo = get(ode_kwargs, :alg, Tsit5()) - sol = solve(ivp, alg=algo; ode_kwargs...) + sol = solve(ivp; alg=algo, ode_kwargs...) T = sol.t - - return reduce(hcat,sol.u) + + return reduce(hcat, sol.u) end if print_times - return jacobian(_flow,backend,λ), T + return jacobian(_flow, backend, λ), T else - return jacobian(_flow,backend,λ) + return jacobian(_flow, backend, λ) end - end - - function ∂λ_flow(t0::Real,x0::Function, tf::Real, λ0::Vector{<:Real}; print_times=false, ode_kwargs...) + end + + function ∂λ_flow( + t0::Real, + x0::Function, + tf::Real, + λ0::Vector{<:Real}; + print_times=false, + ode_kwargs..., + ) T = [] function _flow(λ) - ivp = ODEProblem(rhs, x0(λ), (t0,tf), λ) + ivp = ODEProblem(rhs, x0(λ), (t0, tf), λ) algo = get(ode_kwargs, :alg, Tsit5()) - sol = solve(ivp, alg=algo; ode_kwargs...) + sol = solve(ivp; alg=algo, ode_kwargs...) T = sol.t return sol.u[end] end if print_times - return jacobian(_flow,backend,λ0), T + return jacobian(_flow, backend, λ0), T else - return jacobian(_flow,backend,λ0) + return jacobian(_flow, backend, λ0) end - end - return ∂λ_flow + end + return ∂λ_flow end end # CTDiffFlow diff --git a/src/myode43/figHairer.jl b/src/myode43/figHairer.jl index d699483..253e139 100644 --- a/src/myode43/figHairer.jl +++ b/src/myode43/figHairer.jl @@ -14,53 +14,81 @@ 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₂] +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; +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) +Atol=1.e-4; +Rtol=1.e-4; +x0=[1.5, 3]; +par = [3.0] +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) +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") +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"); +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, +) -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)) +plot(p1, p2, p3; layout=(3, 1)) #= # @@ -90,5 +118,4 @@ 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 index 1ad61db..820eb5a 100644 --- a/src/myode43/myode43.jl +++ b/src/myode43/myode43.jl @@ -4,109 +4,109 @@ # # 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 +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) +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)] # - # 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; + # + # 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 + return T, X +end - function get_Xij(X,i,j) +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] + Xij[l] = X[l][i, j] end return Xij - end - - +end diff --git a/src/myode43/myode43_for_graph.jl b/src/myode43/myode43_for_graph.jl index 11cf063..440ec1f 100644 --- a/src/myode43/myode43_for_graph.jl +++ b/src/myode43/myode43_for_graph.jl @@ -4,106 +4,115 @@ # # 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) +function inith(rhs, x0, par, t0, Atol, Rtol) + # # - # 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; + 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 - 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; + h1 = (0.01/max(d1, d2))^(1/4); + end + h = min(100*h0, h1); + return h +end - return T,X,haccept,Trej,hrej,Ttot,err_loc_est,err_loc_exact,err_abs - 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