1

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:

enter image description here

keepAlive
  • 6,369
  • 5
  • 24
  • 39
B.Gu
  • 21
  • 3
  • What is a 4D-continuous line? An animation of a line moving in a 3d-space ? – keepAlive Jan 12 '19 at 12:46
  • Looks like he wants a 3d line that changes color according to the 4th function in his link. – okovko Jan 12 '19 at 13:17
  • Exactly. I meant continous as opposed to discrete values representations (for which I did find some codes). – B.Gu Jan 13 '19 at 13:09
  • My question wasn't about continuity. What are your dimensions? (Excitation, time, color, ????). What is the 4th? – keepAlive Jan 17 '19 at 20:56
  • @Kanak Thanks for editing the question to make it clearer, I provided a nice solution as a result (: – okovko Jan 18 '19 at 00:29
  • 1
    If any answer does what you want, please consider ticking it as correct. A reputation of 1 is enough to do it. I remind you this because newcomers often forget to do so. See [What should I do when someone answers my question?](https://stackoverflow.com/help/someone-answers) That being said, welcome on [SO](https://stackoverflow.com/). – keepAlive Jan 18 '19 at 16:18
  • Sorry guys, I wasn't available yesterday. This looks great, although the difference between the blue line and the heat (colored) line is still bothering. Definitely enough for me needs though.. Thanks! – B.Gu Jan 19 '19 at 20:24
  • @B.Gu The blue line is what you get from `solve_ivp`. If you know how to get more data points using `solve_ivp`, then go ahead and do so. You could plot the points along the line segments of your original data, if you wanted, but then you wouldn't have a continuous line. – okovko Jan 19 '19 at 20:37

1 Answers1

3

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()

Your curve

okovko
  • 1,851
  • 14
  • 27