How to use pictures in plot to show area preservation using julia?

General Usage
,
1

I have a data of symplectic Euler method

using Plots
h = 0.05;
N = 200;
const g = 9;
const ℓ = 1.46;
θ  = zeros(N+1)
v = zeros(N+1)
θ[1]  = 1.0; # this is θ₀
v[1] = 0.5; # this is v₀
for k in 1:N  
   θ[k+1]  = θ[k] + h * v[k]
   v[k+1]  = v[k] - h * ( g / ℓ ) * sin(θ[k+1])   
 end
plot(θ,v,lc=:red, lw=1)

pend

The flow of this method preserve the area and I would like to investigate this property by using some pictures inside the plot like in the following

Screenshot from 2023-09-06 14-36-36
Screenshot from 2023-09-06 14-37-25

2

Take as starting positions a number of points on a circle and compute the iterates, F^{kn}, of these points, for a fixed n, and k=1, 2, …kmax. Plot the interpolating curve of the corresponding points after each k.

3

Here is the code that implements what is described above:

using Plots
function F(p::Vector{T}; h=0.05, g=9, ℓ=1.46) where T<:AbstractFloat
   Theta =p[1]+h*p[2]
   return [Theta, p[2]-h*(g/ℓ)*sin(Theta)]
end  
composition(F, n) = ∘(ntuple(_ -> F, n)...)

function circle(t; center=[1.0, 0], r=0.3)
    [center[1]+r*cos(t), center[2]+r*sin(t)]
end    

t=range(0, 2pi, 24)
pc=circle.(t;)
pl= plot(getindex.(pc, 1), getindex.(pc, 2), lc=:red, lw=1, 
         legend=false, size=(450, 450), aspect_ratio = :equal)
for _ in 1:10
    pc=composition(F, 5).(pc)
    plot!(pl, getindex.(pc, 1), getindex.(pc, 2), lc=:red, lw=1)
end    
display(pl)
4

Thanks @empet ! The code work perfectly in this case, but when using explicit Euler like

using Plots
h = 0.05;
N = 200;
const g = 9;
const ℓ = 1.46;
θ  = zeros(N+1)
v = zeros(N+1)
θ[1]  = 1.0; # this is θ₀
v[1] = 0.5; # this is v₀
for k in 1:N  
   θ[k+1]  = θ[k] + h * v[k]
   v[k+1]  = v[k] - h * ( g / ℓ ) * sin(θ[k])   
 end
plot(θ,v)

it didn’t give me the expected result like
Screenshot from 2023-09-08 08-19-16

5

Please post your code that generated your last plot.
The function F has the same effect as your discrete system.

7

The first code contain sin(θ[k+1]) and the last code contains sin(θ[k]) which give me
Screenshot from 2023-09-08 10-29-24

8

Without your complete code I cannot help you.

9

Sorry, @empet see this code about explicit Euler

using Plots
function F()
h = 0.05;
N = 200;
 g = 9;
ℓ = 1.46;
θ  = zeros(N+1)
v = zeros(N+1)
θ[1]  = 1.0; # this is θ₀
v[1] = 0.5; # this is v₀
for k in 1:N  
   θ[k+1]  = θ[k] + h * v[k]
   v[k+1]  = v[k] - h * ( g / ℓ ) * sin(θ[k])   
 end
   return hcat(θ, v)
end  

composition(F, n) = ∘(ntuple(_ -> F, n)...)

function circle(t,k)
    [1+(inv(1+h^2))^k*cos(t), (inv(1+h^2))^k*sin(t)]
end    

t=range(0, 2pi, 35)
pc=circle.(t;)
pl= plot(getindex.(pc, 1), getindex.(pc, 2), lc=:red, lw=1, 
         legend=false, size=(450, 450), aspect_ratio = :equal)
for _ in 1:5:100
    pc=composition(F, 10).(pc)
    plot!(pl, getindex.(pc, 1), getindex.(pc, 2), lc=:red, lw=1)
    plot!(x[:,1],x[:,2])
end    
display(pl)
10

Your discrete system defined today is not the same as that you posted initially, when you formulated your question.
The last one, defined today,

θ[k+1]  = θ[k] + h * v[k]
v[k+1]  = v[k] - h * ( g / ℓ ) * sin(θ[k])

is not defined by an area preserving map, because the determinant of its Jacobian matrix is not 1 (one). That’s why your plot reveals an expansive effect, as the number of iterations increases.

Your initial system I have used to illustrate the area preserving of a disk under the action of the map F. was:

 θ[k+1]  = θ[k] + h * v[k]
 v[k+1]  = v[k] - h * ( g / ℓ ) * sin(θ[k+1])
1 Like