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

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.0; # this is θ₀
v = 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)
`````` 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  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.

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+h*p
return [Theta, p-h*(g/ℓ)*sin(Theta)]
end
composition(F, n) = ∘(ntuple(_ -> F, n)...)

function circle(t; center=[1.0, 0], r=0.3)
[center+r*cos(t), center+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)
``````

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.0; # this is θ₀
v = 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 The function F has the same effect as your discrete system.

The first code contain `sin(θ[k+1])` and the last code contains `sin(θ[k])` which give me 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.0; # this is θ₀
v = 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)
``````

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