Plotting a 4-D continous line in Python

Question

I used the solve_ivp function in Pyhon in oredr to solve a system of four coupled differential equations, and now I wish to plot the solution in order to show its stability close to a certain fixed point, in the 4-D phase space.

I was thinking a continous line in a 3-D space, that also changes its color according to the 4th dimension might be a good solution, but couldn't find any way to do so and would love some help.

The code I have used to solve the equations, and plot the four vectors as a function of time:

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

#parametres 
Alpha1= 0.4 #in [0,1]
Alpha2= 1 - Alpha1
Kappa=2 
A1=A2=A=0.3
B1=B2=B=0.03
K1 = 2
K2 = 2
Taun = 3.1
Taup = 3.1

# Differential Equations
def V( t , y):
    dN1dt = Alpha1*Kappa + K1*y[1]*(B1-y[0]) - K2*y[0]*(B2-y[1]) -(y[0]/Taun) - A1*( y[0]*(y[2])-(B1-y[0])*y[2]) 
    dN2dt = Alpha2*Kappa + K2*y[0]*(B2-y[1]) - K1*y[1]*(B1-y[0]) -(y[1]/Taun) - A2*( y[1]*(y[3])-(B2-y[1])*y[3]) 
    dS1dt = -(y[2]/Taup) + A1*(y[0]*(y[2]) - (B1-y[0])*y[2]) 
    dS2dt = -(y[3]/Taup) + A2*(y[1]*(y[3]) - (B2-y[1])*y[3]) 
    return [dN1dt, dN2dt, dS1dt, dS2dt] 



# Stable Points
N1ss = B1/2 + 1/(2*A1*Taup)
N2ss = B2/2 + 1/(2*A2*Taup)
S1ss = Alpha1*Kappa*Taup - 0.5*((B*Taup/Taun) + 1/(A*Taun) + (K1-K2)*Taup*(1/(2*(A**2)*(Taup**2)) - 0.5*(B**2)))
S2ss = Alpha2*Kappa*Taup - 0.5*((B*Taup/Taun) + 1/(A*Taun) + (K2-K1)*Taup*(1/(2*(A**2)*(Taup**2)) - 0.5*(B**2)))
print(N1ss,N2ss,S1ss,S2ss)

# Initial Conditions 
N10 = N1ss*(1.1)
N20 = N2ss*(0.9)
S10 = S1ss*(1.5)
S20 = S2ss*(0.3)
y0 = [N10,N20,S10,S20]

#Solution and Plotting
sol = solve_ivp(V, (0,100) , y0)

plt.plot(sol.t, sol.y.T)
plt.legend(['N1', 'N2', 'S1', 'S2'])
plt.xlabel('Time')
plt.ylabel('Excited Dots and Photons')

Which gives the following plot:

Answer 1

With pyplot, you have to use a scatterplot with many points to draw a colored line. Your data is rather sparse, so you have to fit your parametrized curves and sample them to get more points so you can graph a nice curve. This was fun, I hope this is helpful.

# additional imports

import matplotlib as mpl
from matplotlib import cm as cm
from mpl_toolkits.mplot3d import Axes3D    
from scipy.interpolate import CubicSpline

# include your solving code from above

ts = sol.t
points = sol.y.T
xs = [p[0] for p in points]
ys = [p[1] for p in points]
zs = [p[2] for p in points]
cs = [p[3] for p in points]
colors = cs / max(cs)

csx = CubicSpline(ts, xs)
csy = CubicSpline(ts, ys)
csz = CubicSpline(ts, zs)
csc = CubicSpline(ts, cs)

fine_ts = np.arange(ts[0], ts[-1], 0.01)
fine_xs = [csx(t) for t in fine_ts]
fine_ys = [csy(t) for t in fine_ts]
fine_zs = [csz(t) for t in fine_ts]
fine_cs = [csc(t) for t in fine_ts]
fine_colors = fine_cs / max(fine_cs)

fig = plt.figure(figsize = (10, 10))
ax = fig.gca(projection = '3d')
ax.plot(xs, ys, zs)
ax.scatter(fine_xs, fine_ys, fine_zs, c = cm.hot(fine_colors))
plt.show()