Plot points moving along the velocity field in Julia

Question

The code below creates velocity field. I plotted a blue point in (0.5,0.5). How do I plot series of points that move alongside the velocity field?

using PyPlot

xs = range(0,1,step=0.03)
ys = range(0,1,step=0.03)

nfreq = 20
as = randn(nfreq, nfreq)
aas = randn(nfreq, nfreq)
bs = randn(nfreq, nfreq)
bbs = randn(nfreq, nfreq)

f(x,y) = sum( as[i,j]*sinpi(x*i+ aas[i,j])*sinpi(y*j )/(i^2+j^2)^(1.5) for i=1:nfreq, j=1:nfreq)
g(x,y) = sum( bs[i,j]*sinpi(x*i)*sinpi(y*j + bbs[i,j])/(i^2+j^2)^(1.5) for i=1:nfreq, j=1:nfreq)

quiver(xs,ys, f.(xs,ys'), g.(xs,ys'))

plot

Answer 1

You could use DifferentialEquations.jl, specifically one of the ODE solvers from OrdinaryDiffEq.jl.

Here is a basic (but very useful) tutorial.

In this case, to compute the trajectory starting at [0.1, 0.1], we could do something like:

using OrdinaryDiffEq

F(x, y) = [f(x, y), g(x, y)]
F(u) = F(u...)
# For ODE solver
F(u, t, p) = F(u)

u0 = [0.1, 0.1]
tspan = (0.0, 1.0)
prob = OrdinaryDiffEq.ODEProblem(F, u0, tspan)
sol = OrdinaryDiffEq.solve(prob, OrdinaryDiffEq.Tsit5())

This uses one of the standard solvers to compute the trajectory. The F(u,t,p) method is defined to comply with the signature expected by ODEProblem.

The output sol is a special data structure that stores the (automatically-chosen) time points, the trajectory/solution data, and some info about how the computation went:

julia> sol

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 16-element Vector{Float64}:
 0.0
 0.016059769627848473
 0.05687861311174197
 0.11905519520284545
 0.1762669960378554
 0.24797417284787635
 0.3389099429346001
 0.4168117846405069
 0.5058131414831007
 0.5832877232134522
 0.6673494607419321
 0.7443426334531882
 0.8182342011387053
 0.89201972618978
 0.9834140190395773
 1.0
u: 16-element Vector{Vector{Float64}}:
 [0.1, 0.1]
 [0.10110680137385687, 0.10198760564184167]
 [0.10406194191642619, 0.10713626068506982]
 [0.10896217507542232, 0.11524940132731033]
 [0.11390233972543964, 0.12300643621930929]
 [0.1206820706605841, 0.13312809220811475]
 [0.13023297520058397, 0.14662230992549455]
 [0.1392904751559556, 0.15881027331987904]
 [0.15068267631556032, 0.17352625961041768]
 [0.16157943680595235, 0.18709637923697067]
 [0.17457595399222012, 0.20265945610197006]
 [0.1877136474604169, 0.21764486386096143]
 [0.2015199067922626, 0.23261870623234943]
 [0.21645294590353964, 0.2481557501686445]
 [0.236427795874632, 0.26846545434236724]
 [0.24021585799812886, 0.27231898207853394]

To access the trajectory itself, use sol.u:

julia> sol.u
16-element Vector{Vector{Float64}}:
 [0.1, 0.1]
 [0.10110680137385687, 0.10198760564184167]
 [0.10406194191642619, 0.10713626068506982]
 [0.10896217507542232, 0.11524940132731033]
 [0.11390233972543964, 0.12300643621930929]
 [0.1206820706605841, 0.13312809220811475]
 [0.13023297520058397, 0.14662230992549455]
 [0.1392904751559556, 0.15881027331987904]
 [0.15068267631556032, 0.17352625961041768]
 [0.16157943680595235, 0.18709637923697067]
 [0.17457595399222012, 0.20265945610197006]
 [0.1877136474604169, 0.21764486386096143]
 [0.2015199067922626, 0.23261870623234943]
 [0.21645294590353964, 0.2481557501686445]
 [0.236427795874632, 0.26846545434236724]
 [0.24021585799812886, 0.27231898207853394]

You could then plot this data as a scatter plot along with the vector field itself.