Kernel: SageMath (Development)
Solving a differential equation with SageMath
In [1]:
version()
'SageMath version 8.2.rc2, Release Date: 2018-04-10'
In [2]:
%display latex
Exact solutions: desolve
In [3]:
x = var('x') y = function('y')(x) y
y(x)
In [4]:
eq = diff(y, x) - y == x*y^4 eq
−y(x)+∂x∂y(x)=xy(x)4
In [5]:
print(eq)
-y(x) + diff(y(x), x) == x*y(x)^4
In [6]:
desolve(eq, y)
(−31(3x−1)e(3x)+C)31ex
In [7]:
print(desolve(eq, y))
e^x/(-1/3*(3*x - 1)*e^(3*x) + _C)^(1/3)
In [8]:
desolve(eq, y, show_method='True')
(−31(3x−1)e(3x)+C)31ex,bernoulli
In [9]:
desolve(eq, y, ics=[0, 2])
(−31(3x−1)e(3x)−245)31ex
In [10]:
f(x) = desolve(eq, y, ics=[0, 2]) diff(f(x), x) - f(x) - x*f(x)^4
−(−31(3x−1)e(3x)−245)34xe(4x)+3(−31(3x−1)e(3x)−245)34((3x−1)e(3x)+e(3x))ex
In [11]:
z = diff(f(x), x) - f(x) - x*f(x)^4 z
−(−31(3x−1)e(3x)−245)34xe(4x)+3(−31(3x−1)e(3x)−245)34((3x−1)e(3x)+e(3x))ex
In [12]:
z = diff(f(x), x) - f(x) - x*f(x)^4
In [13]:
z.simplify_full()
0
In [14]:
f(0)
2
System of differential equations
In [15]:
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))
In [16]:
A = matrix([[2,-2,0], [-2,0,2], [0,2,2]]) A
2−20−202022
In [17]:
eqs = [diff(y[i],x) == (A*y)[i] for i in range(3)] eqs
[∂x∂y1(x)=2y1(x)−2y2(x),∂x∂y2(x)=−2y1(x)+2y3(x),∂x∂y3(x)=2y2(x)+2y3(x)]
In [18]:
for eq in eqs: show(eq)
∂x∂y1(x)=2y1(x)−2y2(x)
∂x∂y2(x)=−2y1(x)+2y3(x)
∂x∂y3(x)=2y2(x)+2y3(x)
In [19]:
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)]
Numerical solutions
In [20]:
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')
In [21]:
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)
∂r∂m(r)=4πr2ρ(r)
∂r∂Φ(r)=r2Gm(r)
∂r∂p(r)=−r2Gm(r)ρ(r)
In [22]:
k = var('k') gam = var('gam', latex_name=r'\gamma') p_eos(r) = k*rho(r)^gam p_eos(r)
kρ(r)γ
In [23]:
eq3_rho = eq3.substitute_function(p, p_eos) eq3_rho
γkρ(r)γ−1∂r∂ρ(r)=−r2Gm(r)ρ(r)
In [24]:
eq3_rho = (eq3_rho / (gam*k*rho(r)^(gam-1))).simplify_full() eq3_rho
∂r∂ρ(r)=−γkr2Gρ(r)−γ+2m(r)
In [25]:
eqs = [eq1, eq2, eq3_rho] for eq in eqs: show(eq)
∂r∂m(r)=4πr2ρ(r)
∂r∂Φ(r)=r2Gm(r)
∂r∂ρ(r)=−γkr2Gρ(r)−γ+2m(r)
In [26]:
k0 = 1/4 gam0 = 2 rhs = [eq.rhs().subs({k: k0, gam: gam0, G: 1}) for eq in eqs] rhs
[4πr2ρ(r),r2m(r),−r22m(r)]
In [27]:
rhs[0] = rhs[0] * unit_step(rho(r)) rhs[2] = rhs[2] * unit_step(rho(r)) rhs
[4πr2ρ(r)u(ρ(r)),r2m(r),−r22m(r)u(ρ(r))]
In [28]:
rhs.append(1 * unit_step(rho(r))) rhs
[4πr2ρ(r)u(ρ(r)),r2m(r),−r22m(r)u(ρ(r)),u(ρ(r))]
In [29]:
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),r2m1,−r22m1u(ρ1),u(ρ1)]
In [30]:
rho_c = 1 r_min = 1e-8 r_max = 1 np = 200 delta_r = (r_max - r_min) / (np-1)
In [47]:
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)
In [48]:
delta_r
0.00502512557788945
In [49]:
sol[:10]
[[1.00000000000000×10−8,0,0,1,1.00000000000000×10−8],[0.00502513557789,5.31450779388×10−07,6.61093383901×10−05,0.999867781323,0.00502513557789],[0.0100502611558,4.25098772067×10−06,0.000225597211871,0.999548805576,0.0100502611558],[0.0150753867337,1.43426583858×10−05,0.000490131925791,0.999019736148,0.0150753867337],[0.0201005123116,3.39824144478×10−05,0.000860123055771,0.998279753888,0.0201005123116],[0.0251256378894,6.63340827692×10−05,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]]
In [50]:
rho_sol = [(s[0], s[3]) for s in sol] rho_sol[:10]
[(1.00000000000000×10−8,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)]
In [51]:
graph = line(rho_sol, axes_labels=[r'$r$', r'$\rho$'], gridlines=True) graph
In [52]:
Phi_sol = [(s[0], s[2]) for s in sol] Phi_sol[:10]
[(1.00000000000000×10−8,0),(0.00502513557789,6.61093383901×10−05),(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)]
In [53]:
graph = line(Phi_sol, axes_labels=[r'$r$', r'$\Phi$'], gridlines=True) graph
In [54]:
M = Manifold(2, 'M') print(M)
2-dimensional differentiable manifold M
In [55]:
X.<x,y> = M.chart()
In [56]:
X
(M,(x,y))
In [57]:
v = 2*x*X.frame()[0] - X.frame()[1]
In [58]:
v.display()
2x∂x∂−∂y∂
In [59]:
v.plot()
In [60]:
import numpy as np
In [61]:
z = cos(x).series(x, 8)
In [62]:
z.operands()
[1,−21x2,241x4,−7201x6,O(x8)]
In [0]:
z?
In [0]:
In [0]:
In [0]: