How to know the equation of the plot graphs?

Question

I plotted the 1st-order non-linear differential equations. But I don't know how to know the equations of the plots. Please forgive my ignorance.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import math

def dy_dx(y,x,z):
    c_1 = 5.0 / (1.38 * 1223.0 * pow(10.0, 28.0)*pow(z,3.0))
    c_2 = pow(10.0, 5.0)
    return c_1 * ( y/math.sqrt(1.0+pow(y, 2.0)) ) * ( ((1.0-pow(y, 3.0))/(z* pow(y,(1.0/3.0)))) - (y * c_2) )

xs = np.linspace(0, pow(10.0, 12.0), pow(10.0, 6.0))
y_0 = 1.0
z = 0.00001
y1 = odeint(dy_dx, y_0, xs, args=(z,))
z = 0.000015
y2 = odeint(dy_dx, y_0, xs, args=(z,))
y1 = np.array(y1).flatten()
y2 = np.array(y2).flatten()

plt.rcParams.update({'font.size': 10})
plt.ylim(0,1.0)
plt.xlabel("x")
plt.ylabel("y")
plt.plot(xs, y1, 'r-')
plt.plot(xs, y2, 'b-')
plt.grid(True)
plt.show()

Answer 1

It is very improbable that any symbolic solution exists, with third powers and third roots on the right side. The relevant factors for the qualitative behavior are

y^(2/3)*(1-y^3-c2*z*y^(4/3))

what remains are positive non-zero factors. The one-dimensional dynamic is very simple, roots of the right side are equilibrium positions, the sign in-between tells you if the solutions in that interval are falling or growing.

• At y=0 you can have branching of solutions similar to y'=C*y^(2/3), as the function is not locally Lipschitz there.

• The polynomial u^9+a*u^4-1=0, a=c_2*z has exactly one positive root u_ast, y_equi=u_ast^3, for small a it is close to u=1 (next approximation u=1/(1+a)^(1/4)~1/(1+a/4)), for large a close to u=a^(-1/4), y=a^(-3/4); and zero or two negative roots. As the right side of the ODE changes sign from positive to negative when passing over y_equi, this is a stable equilibrium.

  In your examples with c2=1e5 and z=1e-5 or z=1.5e-5 one gets a=1 or a=1.5 which is in-between, the intermediate value theorem tells us that there is a root between 0 and 1, u_ast=0.89345331 and y_equi=0.71320698 in the first case and u_ast=0.84759991 and y_equi=0.60893748 in the second. This is exactly what you see in the plot.

That is about all of the theoretical content that you can hope to extract from this equation for positive y.