Phase portrait with vector field

[moxeeeem]

I’m trying to reformat a python code that plots a phase portrait with vector field. It seems that I did everything right, but the error is in the quiver () function. What did I write wrong?

Python code (working correctly):

# Phase portrait with vector field. Check two systems are the same.
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
import pylab as pl

# The 2-dimensional linear system.
a, b, c, d = 2, 1, 1, 2
def dx_dt(x, t):
    return [a*x[0] + b*x[1], c*x[0] + d*x[1]]

# Trajectories in forward time.
ts = np.linspace(0, 4, 100)
ic = np.linspace(-1, 1, 5)
for r in ic:
    for s in ic:
        x0 = [r, s]
        xs = odeint(dx_dt, x0, ts)
        plt.plot(xs[:,0], xs[:,1], "r-")

# Trajectories in backward time.
ts = np.linspace(0, -4, 100)
ic = np.linspace(-1, 1, 5)
for r in ic:
    for s in ic:
        x0 = [r, s]
        xs = odeint(dx_dt, x0, ts)
        plt.plot(xs[:,0], xs[:,1], "r-")
        
# Label the axes and set fontsizes.
plt.xlabel("x", fontsize=15)
plt.ylabel("y", fontsize=15)
plt.tick_params(labelsize=15)
plt.xlim(-1, 1)
plt.ylim(-1, 1);

# Plot the vectorfield.

X,Y = np.mgrid[-1:1:10j, -1:1:10j]
u = a*X + b*Y
v = c*X + d*Y
pl.quiver(X, Y, u, v, color = "b")
plt.show()

Output which I want to get:

My Julia code:

using Plots
using SciPy
using LinearAlgebra
using DifferentialEquations

# The 2-dimensional linear system.
a, b, c, d = 2, 1, 1, 2

function dx_dt(du, u, p, t)
    du[1] = a*u[1] + b*u[2]
    du[2] = c*u[1] + d*u[2]
end

# Trajectories in forward time.
ts = LinRange(0, 4, 100)
ic = LinRange(-1, 1, 5)
tspan = (0.0, 4.0)

for r in ic
    for s in ic
        x0 = [r, s]
        xs = solve(ODEProblem(dx_dt, x0, tspan))
        Plots.plot(xs, color="red", linestyle="-")
        
    end
end

# Trajectories in backward time.
ts = LinRange(0, -4, 100)
ic = LinRange(-1, 1, 5)
tspan = (-4.0, 0.0)  

for r in ic
    for s in ic
        x0 = [r, s]
        xs = solve(ODEProblem(dx_dt, x0, tspan))
        Plots.plot(xs, color="red", linestyle="-")
    end
end
        
# Label the axes and set fontsizes.
xlabel!("x")
ylabel!("y")
#fontfamily("sans-serif")
#fontsize!(15)
xlims!(-1, 1)
ylims!(-1, 1)

# Plot the vectorfield.
function meshgrid(x, y)
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    return (X, Y)
end

x = range(-1, 1, length=10)
y = range(-1, 1, length=10)

a, b, c, d = 1, 1, -1, 1

X, Y = meshgrid(x, y)

u = a .* X .+ b .* Y
v = c .* X .+ d .* Y

quiver(X, Y, u, v, color="blue", alpha=0.5, aspect_ratio=:equal, framestyle=:box)```

My error: 

Couldn’t process recipe args: (Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64})

[tomerarnon]

It’s usually a mistake to make a “meshgrid” in julia…

quiver takes in 4 vectors, and has signature quiver(x, y, quiver = (u, v)). x and y are vectors of x-y coordinates. It will also accept a vector of 2-Tuples (each one being a point), or an Nx2 matrix. u and v must also be vectors, representing the x and y offset of each arrow in the quiver. You have all 4 being NxN matrices.

Given how you’ve set things up (with a meshgrid), you could recover the vectors you need by calling vec on x, y, u and v. In general though, it is better not to use a meshgrid in the first place and just iterate over the product of x and y with e.g. iterators.product(x, y), or using a nested for loop, or something similar.

Using vec, I get what I assume you were looking for:
quiver(vec(X), vec(Y), quiver=(vec(u), vec(v)), color="blue", alpha=0.5, aspect_ratio=:equal, framestyle=:box)

[moxeeeem]

Thank you very much, your answer saved me. Do you know why I don’t have the red lines of my plot displayed?

[tomerarnon]

Not off the top of my head, because I haven’t gone over all of your code in detail. What I did notice is this:

• you use Plots.plot repeatedly. Each of these makes a new plot. To add things to an existing plot, you need to use the plot! form (or quiver!, etc)
• I don’t think linestyle="-" is something thats supported… not sure where you got that (I could be wrong about this though)

I suggest splitting this up and solving this one piece at a time. Make a plot with the red lines separately from the quiver, and combine once you figure out each piece.