Shareddiff_solve_tour.ipynbOpen in CoCalc

Solving a differential equation with SageMath

version()

'SageMath version 8.2.rc2, Release Date: 2018-04-10'
%display latex


Exact solutions: desolve

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

$y\left(x\right)$
eq = diff(y, x) - y == x*y^4
eq

$-y\left(x\right) + \frac{\partial}{\partial x}y\left(x\right) = x y\left(x\right)^{4}$
print(eq)

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

$\frac{e^{x}}{{\left(-\frac{1}{3} \, {\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} + C\right)}^{\frac{1}{3}}}$
print(desolve(eq, y))

e^x/(-1/3*(3*x - 1)*e^(3*x) + _C)^(1/3)
desolve(eq, y, show_method='True')

$\left[\frac{e^{x}}{{\left(-\frac{1}{3} \, {\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} + C\right)}^{\frac{1}{3}}}, \verb|bernoulli|\right]$
desolve(eq, y, ics=[0, 2])

$\frac{e^{x}}{{\left(-\frac{1}{3} \, {\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} - \frac{5}{24}\right)}^{\frac{1}{3}}}$
f(x) = desolve(eq, y, ics=[0, 2])
diff(f(x), x) - f(x) - x*f(x)^4

$-\frac{x e^{\left(4 \, x\right)}}{{\left(-\frac{1}{3} \, {\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} - \frac{5}{24}\right)}^{\frac{4}{3}}} + \frac{{\left({\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} + e^{\left(3 \, x\right)}\right)} e^{x}}{3 \, {\left(-\frac{1}{3} \, {\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} - \frac{5}{24}\right)}^{\frac{4}{3}}}$
z = diff(f(x), x) - f(x) - x*f(x)^4
z

$-\frac{x e^{\left(4 \, x\right)}}{{\left(-\frac{1}{3} \, {\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} - \frac{5}{24}\right)}^{\frac{4}{3}}} + \frac{{\left({\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} + e^{\left(3 \, x\right)}\right)} e^{x}}{3 \, {\left(-\frac{1}{3} \, {\left(3 \, x - 1\right)} e^{\left(3 \, x\right)} - \frac{5}{24}\right)}^{\frac{4}{3}}}$
z = diff(f(x), x) - f(x) - x*f(x)^4

z.simplify_full()

$0$
f(0)

$2$

System of differential equations

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

$\left(y_{1}\left(x\right),\,y_{2}\left(x\right),\,y_{3}\left(x\right)\right)$
A = matrix([[2,-2,0], [-2,0,2], [0,2,2]])
A

$\left(\begin{array}{rrr} 2 & -2 & 0 \\ -2 & 0 & 2 \\ 0 & 2 & 2 \end{array}\right)$
eqs = [diff(y[i],x) == (A*y)[i] for i in range(3)]
eqs

$\left[\frac{\partial}{\partial x}y_{1}\left(x\right) = 2 \, y_{1}\left(x\right) - 2 \, y_{2}\left(x\right), \frac{\partial}{\partial x}y_{2}\left(x\right) = -2 \, y_{1}\left(x\right) + 2 \, y_{3}\left(x\right), \frac{\partial}{\partial x}y_{3}\left(x\right) = 2 \, y_{2}\left(x\right) + 2 \, y_{3}\left(x\right)\right]$
for eq in eqs:
show(eq)

$\frac{\partial}{\partial x}y_{1}\left(x\right) = 2 \, y_{1}\left(x\right) - 2 \, y_{2}\left(x\right)$
$\frac{\partial}{\partial x}y_{2}\left(x\right) = -2 \, y_{1}\left(x\right) + 2 \, y_{3}\left(x\right)$
$\frac{\partial}{\partial x}y_{3}\left(x\right) = 2 \, y_{2}\left(x\right) + 2 \, y_{3}\left(x\right)$
sol = desolve_system(eqs, [y1,y2,y3], ics=[0, 2, 1, -2])
sol

$\left[y_{1}\left(x\right) = e^{\left(4 \, x\right)} + e^{\left(-2 \, x\right)}, y_{2}\left(x\right) = -e^{\left(4 \, x\right)} + 2 \, e^{\left(-2 \, x\right)}, y_{3}\left(x\right) = -e^{\left(4 \, x\right)} - e^{\left(-2 \, x\right)}\right]$

Numerical solutions

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)

$\frac{\partial}{\partial r}m\left(r\right) = 4 \, \pi r^{2} \rho\left(r\right)$
$\frac{\partial}{\partial r}\Phi\left(r\right) = \frac{G m\left(r\right)}{r^{2}}$
$\frac{\partial}{\partial r}p\left(r\right) = -\frac{G m\left(r\right) \rho\left(r\right)}{r^{2}}$
k = var('k')
gam = var('gam', latex_name=r'\gamma')
p_eos(r) = k*rho(r)^gam
p_eos(r)

$k \rho\left(r\right)^{{\gamma}}$
eq3_rho = eq3.substitute_function(p, p_eos)
eq3_rho

${\gamma} k \rho\left(r\right)^{{\gamma} - 1} \frac{\partial}{\partial r}\rho\left(r\right) = -\frac{G m\left(r\right) \rho\left(r\right)}{r^{2}}$
eq3_rho = (eq3_rho / (gam*k*rho(r)^(gam-1))).simplify_full()
eq3_rho

$\frac{\partial}{\partial r}\rho\left(r\right) = -\frac{G \rho\left(r\right)^{-{\gamma} + 2} m\left(r\right)}{{\gamma} k r^{2}}$
eqs = [eq1, eq2, eq3_rho]
for eq in eqs:
show(eq)

$\frac{\partial}{\partial r}m\left(r\right) = 4 \, \pi r^{2} \rho\left(r\right)$
$\frac{\partial}{\partial r}\Phi\left(r\right) = \frac{G m\left(r\right)}{r^{2}}$
$\frac{\partial}{\partial r}\rho\left(r\right) = -\frac{G \rho\left(r\right)^{-{\gamma} + 2} m\left(r\right)}{{\gamma} k r^{2}}$
k0 = 1/4
gam0 = 2
rhs = [eq.rhs().subs({k: k0, gam: gam0, G: 1}) for eq in eqs]
rhs

$\left[4 \, \pi r^{2} \rho\left(r\right), \frac{m\left(r\right)}{r^{2}}, -\frac{2 \, m\left(r\right)}{r^{2}}\right]$
rhs[0] = rhs[0] * unit_step(rho(r))
rhs[2] = rhs[2] * unit_step(rho(r))
rhs

$\left[4 \, \pi r^{2} \rho\left(r\right) \mathrm{u}\left(\rho\left(r\right)\right), \frac{m\left(r\right)}{r^{2}}, -\frac{2 \, m\left(r\right) \mathrm{u}\left(\rho\left(r\right)\right)}{r^{2}}\right]$
rhs.append(1 * unit_step(rho(r)))
rhs

$\left[4 \, \pi r^{2} \rho\left(r\right) \mathrm{u}\left(\rho\left(r\right)\right), \frac{m\left(r\right)}{r^{2}}, -\frac{2 \, m\left(r\right) \mathrm{u}\left(\rho\left(r\right)\right)}{r^{2}}, \mathrm{u}\left(\rho\left(r\right)\right)\right]$
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

$\left[4 \, \pi r^{2} \rho_{1} \mathrm{u}\left(\rho_{1}\right), \frac{m_{1}}{r^{2}}, -\frac{2 \, m_{1} \mathrm{u}\left(\rho_{1}\right)}{r^{2}}, \mathrm{u}\left(\rho_{1}\right)\right]$
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

$0.00502512557788945$
sol[:10]

$\left[\left[1.00000000000000 \times 10^{-8}, 0, 0, 1, 1.00000000000000 \times 10^{-8}\right], \left[0.00502513557789, 5.31450779388 \times 10^{-07}, 6.61093383901 \times 10^{-05}, 0.999867781323, 0.00502513557789\right], \left[0.0100502611558, 4.25098772067 \times 10^{-06}, 0.000225597211871, 0.999548805576, 0.0100502611558\right], \left[0.0150753867337, 1.43426583858 \times 10^{-05}, 0.000490131925791, 0.999019736148, 0.0150753867337\right], \left[0.0201005123116, 3.39824144478 \times 10^{-05}, 0.000860123055771, 0.998279753888, 0.0201005123116\right], \left[0.0251256378894, 6.63340827692 \times 10^{-05}, 0.00133551784062, 0.997328964319, 0.0251256378894\right], \left[0.0301507634673, 0.000114545342118, 0.00191616015049, 0.996167679699, 0.0301507634673\right], \left[0.0351758890452, 0.000181743715046, 0.00260183812477, 0.99479632375, 0.0351758890452\right], \left[0.0402010146231, 0.000271032578936, 0.0033922948088, 0.993215410382, 0.0402010146231\right], \left[0.045226140201, 0.000385487200038, 0.00428723139284, 0.991425537214, 0.045226140201\right]\right]$
rho_sol = [(s[0], s[3]) for s in sol]
rho_sol[:10]

$\left[\left(1.00000000000000 \times 10^{-8}, 1\right), \left(0.00502513557789, 0.999867781323\right), \left(0.0100502611558, 0.999548805576\right), \left(0.0150753867337, 0.999019736148\right), \left(0.0201005123116, 0.998279753888\right), \left(0.0251256378894, 0.997328964319\right), \left(0.0301507634673, 0.996167679699\right), \left(0.0351758890452, 0.99479632375\right), \left(0.0402010146231, 0.993215410382\right), \left(0.045226140201, 0.991425537214\right)\right]$
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]

$\left[\left(1.00000000000000 \times 10^{-8}, 0\right), \left(0.00502513557789, 6.61093383901 \times 10^{-05}\right), \left(0.0100502611558, 0.000225597211871\right), \left(0.0150753867337, 0.000490131925791\right), \left(0.0201005123116, 0.000860123055771\right), \left(0.0251256378894, 0.00133551784062\right), \left(0.0301507634673, 0.00191616015049\right), \left(0.0351758890452, 0.00260183812477\right), \left(0.0402010146231, 0.0033922948088\right), \left(0.045226140201, 0.00428723139284\right)\right]$
graph = line(Phi_sol, axes_labels=[r'$r$', r'$\Phi$'], gridlines=True)
graph

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

2-dimensional differentiable manifold M
X.<x,y> = M.chart()

X

$\left(M,(x, y)\right)$
v = 2*x*X.frame()[0] - X.frame()[1]

v.display()

$2 \, x \frac{\partial}{\partial x } -\frac{\partial}{\partial y }$
v.plot()

import numpy as np

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

z.operands()

$\left[1, -\frac{1}{2} \, x^{2}, \frac{1}{24} \, x^{4}, -\frac{1}{720} \, x^{6}, \mathcal{O}\left(x^{8}\right)\right]$
z?