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(x)y\left(x\right)
eq = diff(y, x) - y == x*y^4
eq
y(x)+xy(x)=xy(x)4-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)
ex(13(3x1)e(3x)+C)13\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')
[ex(13(3x1)e(3x)+C)13,bernoulli]\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])
ex(13(3x1)e(3x)524)13\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
xe(4x)(13(3x1)e(3x)524)43+((3x1)e(3x)+e(3x))ex3(13(3x1)e(3x)524)43-\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
xe(4x)(13(3x1)e(3x)524)43+((3x1)e(3x)+e(3x))ex3(13(3x1)e(3x)524)43-\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()
00
f(0)
22

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
(y1(x),y2(x),y3(x))\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
(220202022)\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
[xy1(x)=2y1(x)2y2(x),xy2(x)=2y1(x)+2y3(x),xy3(x)=2y2(x)+2y3(x)]\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)
xy1(x)=2y1(x)2y2(x)\frac{\partial}{\partial x}y_{1}\left(x\right) = 2 \, y_{1}\left(x\right) - 2 \, y_{2}\left(x\right)
xy2(x)=2y1(x)+2y3(x)\frac{\partial}{\partial x}y_{2}\left(x\right) = -2 \, y_{1}\left(x\right) + 2 \, y_{3}\left(x\right)
xy3(x)=2y2(x)+2y3(x)\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
[y1(x)=e(4x)+e(2x),y2(x)=e(4x)+2e(2x),y3(x)=e(4x)e(2x)]\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)
rm(r)=4πr2ρ(r)\frac{\partial}{\partial r}m\left(r\right) = 4 \, \pi r^{2} \rho\left(r\right)
rΦ(r)=Gm(r)r2\frac{\partial}{\partial r}\Phi\left(r\right) = \frac{G m\left(r\right)}{r^{2}}
rp(r)=Gm(r)ρ(r)r2\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ρ(r)γk \rho\left(r\right)^{{\gamma}}
eq3_rho = eq3.substitute_function(p, p_eos)
eq3_rho
γkρ(r)γ1rρ(r)=Gm(r)ρ(r)r2{\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
rρ(r)=Gρ(r)γ+2m(r)γkr2\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)
rm(r)=4πr2ρ(r)\frac{\partial}{\partial r}m\left(r\right) = 4 \, \pi r^{2} \rho\left(r\right)
rΦ(r)=Gm(r)r2\frac{\partial}{\partial r}\Phi\left(r\right) = \frac{G m\left(r\right)}{r^{2}}
rρ(r)=Gρ(r)γ+2m(r)γkr2\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
[4πr2ρ(r),m(r)r2,2m(r)r2]\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
[4πr2ρ(r)u(ρ(r)),m(r)r2,2m(r)u(ρ(r))r2]\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
[4πr2ρ(r)u(ρ(r)),m(r)r2,2m(r)u(ρ(r))r2,u(ρ(r))]\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
[4πr2ρ1u(ρ1),m1r2,2m1u(ρ1)r2,u(ρ1)]\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.005025125577889450.00502512557788945
sol[:10]
[[1.00000000000000×108,0,0,1,1.00000000000000×108],[0.00502513557789,5.31450779388×1007,6.61093383901×1005,0.999867781323,0.00502513557789],[0.0100502611558,4.25098772067×1006,0.000225597211871,0.999548805576,0.0100502611558],[0.0150753867337,1.43426583858×1005,0.000490131925791,0.999019736148,0.0150753867337],[0.0201005123116,3.39824144478×1005,0.000860123055771,0.998279753888,0.0201005123116],[0.0251256378894,6.63340827692×1005,0.00133551784062,0.997328964319,0.0251256378894],[0.0301507634673,0.000114545342118,0.00191616015049,0.996167679699,0.0301507634673],[0.0351758890452,0.000181743715046,0.00260183812477,0.99479632375,0.0351758890452],[0.0402010146231,0.000271032578936,0.0033922948088,0.993215410382,0.0402010146231],[0.045226140201,0.000385487200038,0.00428723139284,0.991425537214,0.045226140201]]\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]
[(1.00000000000000×108,1),(0.00502513557789,0.999867781323),(0.0100502611558,0.999548805576),(0.0150753867337,0.999019736148),(0.0201005123116,0.998279753888),(0.0251256378894,0.997328964319),(0.0301507634673,0.996167679699),(0.0351758890452,0.99479632375),(0.0402010146231,0.993215410382),(0.045226140201,0.991425537214)]\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]
[(1.00000000000000×108,0),(0.00502513557789,6.61093383901×1005),(0.0100502611558,0.000225597211871),(0.0150753867337,0.000490131925791),(0.0201005123116,0.000860123055771),(0.0251256378894,0.00133551784062),(0.0301507634673,0.00191616015049),(0.0351758890452,0.00260183812477),(0.0402010146231,0.0033922948088),(0.045226140201,0.00428723139284)]\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
(M,(x,y))\left(M,(x, y)\right)
v = 2*x*X.frame()[0] - X.frame()[1]
v.display()
2xxy2 \, x \frac{\partial}{\partial x } -\frac{\partial}{\partial y }
v.plot()
import numpy as np
z = cos(x).series(x, 8)
z.operands()
[1,12x2,124x4,1720x6,O(x8)]\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?