Solving a differential equation with SageMath

version()

%display latex

x = var('x')
y = function('y')(x)
y

eq = diff(y, x) - y == x*y^4
eq

print(eq)

desolve(eq, y)

print(desolve(eq, y))

desolve(eq, y, show_method='True')

desolve(eq, y, ics=[0, 2])

f(x) = desolve(eq, y, ics=[0, 2])
diff(f(x), x) - f(x) - x*f(x)^4

z = diff(f(x), x) - f(x) - x*f(x)^4
z

z = diff(f(x), x) - f(x) - x*f(x)^4

z.simplify_full()

f(0)

y1 = function('y_1')(x)
y2 = function('y_2')(x)
y3 = function('y_3')(x)
y = vector([y1, y2, y3])
y

A = matrix([[2,-2,0], [-2,0,2], [0,2,2]])
A

eqs = [diff(y[i],x) == (A*y)[i] for i in range(3)]
eqs

for eq in eqs:
    show(eq)

sol = desolve_system(eqs, [y1,y2,y3], ics=[0, 2, 1, -2])
sol

rho = function('rho', latex_name=r'\rho')
p = function('p')
m = function('m')
Phi = function('Phi', latex_name=r'\Phi')
r = var('r')
G = var('G')

eq1 = diff(m(r), r) == 4*pi*r^2*rho(r)
eq2 = diff(Phi(r), r) == G*m(r)/r^2
eq3 = diff(p(r), r) == -rho(r)*G*m(r)/r^2
for eq in [eq1, eq2, eq3]:
    show(eq)

k = var('k')
gam = var('gam', latex_name=r'\gamma')
p_eos(r) = k*rho(r)^gam
p_eos(r)

eq3_rho = eq3.substitute_function(p, p_eos)
eq3_rho

eq3_rho = (eq3_rho / (gam*k*rho(r)^(gam-1))).simplify_full()
eq3_rho

eqs = [eq1, eq2, eq3_rho]
for eq in eqs:
    show(eq)

k0 = 1/4
gam0 = 2
rhs = [eq.rhs().subs({k: k0, gam: gam0, G: 1}) for eq in eqs]
rhs

rhs[0] = rhs[0] * unit_step(rho(r))
rhs[2] = rhs[2] * unit_step(rho(r))
rhs

rhs.append(1 * unit_step(rho(r)))
rhs

var('m_1 Phi_1 rho_1 r_1')
rhs = [y.subs({m(r): m_1, Phi(r): Phi_1, rho(r): rho_1}) for y in rhs]
rhs

rho_c = 1
r_min = 1e-8
r_max = 1
np = 200
delta_r = (r_max - r_min) / (np-1)

sol = desolve_system_rk4(rhs, vars=(m_1, Phi_1, rho_1, r_1), ivar=r, 
                         ics=[r_min, 0, 0, rho_c, r_min], 
                         end_points=r_max, step=delta_r)

delta_r

sol[:10]

rho_sol = [(s[0], s[3]) for s in sol]
rho_sol[:10]

graph = line(rho_sol, axes_labels=[r'$r$', r'$\rho$'], gridlines=True)
graph

Phi_sol = [(s[0], s[2]) for s in sol]
Phi_sol[:10]

graph = line(Phi_sol, axes_labels=[r'$r$', r'$\Phi$'], gridlines=True)
graph

M = Manifold(2, 'M')
print(M)

X.<x,y> = M.chart()

X

v = 2*x*X.frame()[0] - X.frame()[1]

v.display()

v.plot()

import numpy as np

z = cos(x).series(x, 8)

z.operands()

z?