Ask Your Question
1

Differential equation solution

asked 2024-02-20 18:19:47 +0100

roux_de_secours gravatar image

I'm still a bit new to sage math and I don't understand the result given by desolve_system, here is my code

t = var('t')
var('w C1 C2 C3 n ')
H = function('H')(t)
phi = (C2 * t + C3)^n #function('\\phi')(t)
sigma = function('\\Sigma')(t)
V = C*phi^2#function('V')(phi)
P = function('P')(t)

######

de1 = H^2-w/6 * (P/phi)^2 - V/(6*phi) - 2*sigma/3 + H*P/phi
de2 = diff(H, t) + w/6 * (P/phi)^2 - V/(3*phi) + 2*sigma/3 + 2*H^2 - 1/(2*phi) * (phi * diff(V, t) - 2*V)/(2*w+3)
de3 = diff(phi, t) == P

des = [de1, de2, de3]

######

vars = [H, sigma, P]
sol= desolve_system(des, vars, ivar = t)

And it gives solutions of the form

$$\newcommand{\Bold}[1]{\mathbf{#1}}P\left(t\right) = \mathcal{L}^{-1}\left(\frac{C_{2}^{n} n e^{\left(\frac{C_{3} g_{63262}}{C_{2}}\right)} \Gamma\left(n, \frac{C_{3} g_{63262}}{C_{2}}\right)}{g_{63262}^{n}}, g_{63262}, t\right)$$

I'm not so sure what is the $\mathcal{L}$ in this case, nor the $g_{63262}$ terms. What am I looking at. I have found no information online, maybe because I wasn't sure how to describe this properly.

Thank you.

edit retag flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted
1

answered 2024-02-21 11:25:04 +0100

Emmanuel Charpentier gravatar image

updated 2024-02-21 14:09:26 +0100

After running your code, typing your solution in text mode gives :

sage: sol
[H(t) == ilt(1/2*(2*(2*w*(H(0) - 3*laplace(H(t)^2, t, g1626) - laplace(H(t)*P(t)/(C2*t + C3)^n, t, g1626)) + 3*H(0) - 9*laplace(H(t)^2, t, g1626) - 3*laplace(H(t)*P(t)/(C2*t + C3)^n, t, g1626))*g1626^(2*n + 1) + (2*C*C2^(2*n)*g1626*n*gamma(2*n, C3*g1626/C2) + (2*C*C2^n*w*gamma(n + 1, C3*g1626/C2) + C*C2^n*gamma(n + 1, C3*g1626/C2))*g1626^n)*e^(C3*g1626/C2))*g1626^(-2*n - 2)/(2*w + 3), g1626, t),
 sigma(t) == ilt(-1/4*(C*C2^n*e^(C3*g1626/C2)*gamma(n + 1, C3*g1626/C2) + (w*laplace(P(t)^2/(C2*t + C3)^(2*n), t, g1626) - 6*laplace(H(t)^2, t, g1626) - 6*laplace(H(t)*P(t)/(C2*t + C3)^n, t, g1626))*g1626^(n + 1))*g1626^(-n - 1), g1626, t),
 P(t) == ilt(C2^n*n*e^(C3*g1626/C2)*gamma(n, C3*g1626/C2)/g1626^n, g1626, t)]

which should lead you to :

sage: sol[0].operands()[1].operator()
ilt
sage: ilt
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[23], line 1
----> 1 ilt

NameError: name 'ilt' is not defined

This is a (long-standing) implementation defect. ilt is Maxima's inverse Laplace transform operator ; from maxima_calculus.ilt? :

-- Function: ilt (<expr>, <s>, <t>)

     Computes the inverse Laplace transform of <expr> with respect
     to <s> and parameter <t>.  <expr> must be a ratio of
     polynomials whose denominator has only linear and quadratic
     factors; there is an extension of ‘ilt’, called ‘pwilt’
     (Piece-Wise Inverse Laplace Transform) that handles several
     other cases where ‘ilt’ fails.

The gxxxx symbols, generated by Maxima's solver, denote integration constants, to be solved to satisfy boundary conditions.

By the way, the first and third equations of your system are ordinary equations ; the only differential equation of your system is the second :

$$ H\left(t\right)^{2} - \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{6 \, {\left(C_{2} t + C_{3}\right)}^{n}} + \frac{H\left(t\right) P\left(t\right)}{{\left(C_{2} t + C_{3}\right)}^{n}} - \frac{w P\left(t\right)^{2}}{6 \, {\left(C_{2} t + C_{3}\right)}^{2 \, n}} - \frac{2}{3} \, \Sigma\left(t\right) = 0 $$

$$ 2 \, H\left(t\right)^{2} - \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{3 \, {\left(C_{2} t + C_{3}\right)}^{n}} + \frac{w P\left(t\right)^{2}}{6 \, {\left(C_{2} t + C_{3}\right)}^{2 \, n}} - \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n - 1} {\left(C_{2} t + C_{3}\right)}^{n} C C_{2} n - {\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{{\left(C_{2} t + C_{3}\right)}^{n} {\left(2 \, w + 3\right)}} + \frac{2}{3} \, \Sigma\left(t\right) + \frac{\partial}{\partial t}H\left(t\right) = 0 $$

$$ {\left(C_{2} t + C_{3}\right)}^{n - 1} C_{2} n = P\left(t\right) $$

We can solve the ordinary equations for $P(t)$ an $\Sigma(t)$ (by temporarily substituting variable for these expressions) and obtain a form f the differential equation where they no longer appear :

 D1=dict(zip((H(t=t), sigma(t=t), P(t=t)), var("f", n=3)))
Sol1=solve([des[u].subs(D1) for u in (0, 2)], (f1, f2))
ID1={D1[u]:u for u in D1.keys()}
DE=des[1].subs(D1).subs(Sol1[0]).subs(ID1)

$$ \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n - 2} C_{2}^{2} n^{2} w}{6 \, {\left(C_{2} t + C_{3}\right)}^{2 \, n}} + 2 \, H\left(t\right)^{2} - \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{3 \, {\left(C_{2} t + C_{3}\right)}^{n}} + \frac{6 \, C_{2}^{2} t^{2} H\left(t\right)^{2} - C_{2}^{2} n^{2} w + 6 \, C_{2} C_{3} n H\left(t\right) + 6 \, C_{3}^{2} H\left(t\right)^{2} - {\left(C C_{2}^{2} t^{2} + 2 \, C C_{2} C_{3} t + C C_{3}^{2}\right)} {\left(C_{2} t + C_{3}\right)}^{n} + 6 \, {\left(C_{2}^{2} n H\left(t\right) + 2 \, C_{2} C_{3} H\left(t\right)^{2}\right)} t}{6 \, {\left(C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}\right)}} - \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n - 1} {\left(C_{2} t + C_{3}\right)}^{n} C C_{2} n - {\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{{\left(C_{2} t + C_{3}\right)}^{n} {\left(2 \, w + 3\right)}} + \frac{\partial}{\partial t}H\left(t\right) = 0 $$

We can rewrite this differential equation by isolating the differential term, which, after resubstituting to a variable is equal to :

DE1=((-DE)+diff(H(t=t), t)).lhs().subs(D1)

$$ -\frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n - 2} C_{2}^{2} n^{2} w}{6 \, {\left(C_{2} t + C_{3}\right)}^{2 \, n}} - 2 \, f_{0}^{2} + \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{3 \, {\left(C_{2} t + C_{3}\right)}^{n}} - \frac{6 \, C_{2}^{2} f_{0}^{2} t^{2} - C_{2}^{2} n^{2} w + 6 \, C_{3}^{2} f_{0}^{2} + 6 \, C_{2} C_{3} f_{0} n - {\left(C C_{2}^{2} t^{2} + 2 \, C C_{2} C_{3} t + C C_{3}^{2}\right)} {\left(C_{2} t + C_{3}\right)}^{n} + 6 \, {\left(2 \, C_{2} C_{3} f_{0}^{2} + C_{2}^{2} f_{0} n\right)} t}{6 \, {\left(C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}\right)}} + \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n - 1} {\left(C_{2} t + C_{3}\right)}^{n} C C_{2} n - {\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{{\left(C_{2} t + C_{3}\right)}^{n} {\left(2 \, w + 3\right)}} $$

This equation is polynomial in $f0$ ; we can rewrite it as :

DE2=diff(H(t=t), t)==sum([u[0]*H(t=t)^u[1] for u in DE1.coefficients(f0)])

$$ \frac{\partial}{\partial t}H\left(t\right) = \frac{{\left(C_{2} t + C_{3}\right)}^{n} C C_{2}^{2} t^{2}}{6 \, {\left(C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}\right)}} + \frac{{\left(C_{2} t + C_{3}\right)}^{n} C C_{2} C_{3} t}{3 \, {\left(C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}\right)}} + \frac{C_{2}^{2} n^{2} w}{6 \, {\left(C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}\right)}} - \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n - 2} C_{2}^{2} n^{2} w}{6 \, {\left(C_{2} t + C_{3}\right)}^{2 \, n}} + \frac{{\left(C_{2} t + C_{3}\right)}^{n} C C_{3}^{2}}{6 \, {\left(C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}\right)}} - {\left(\frac{C_{2}^{2} t^{2}}{C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}} + \frac{2 \, C_{2} C_{3} t}{C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}} + \frac{C_{3}^{2}}{C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}} + 2\right)} H\left(t\right)^{2} - {\left(\frac{C_{2}^{2} n t}{C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}} + \frac{C_{2} C_{3} n}{C_{2}^{2} t^{2} + 2 \, C_{2} C_{3} t + C_{3}^{2}}\right)} H\left(t\right) + \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{3 \, {\left(C_{2} t + C_{3}\right)}^{n}} + \frac{{\left(C_{2} t + C_{3}\right)}^{2 \, n - 1} {\left(C_{2} t + C_{3}\right)}^{n} C C_{2} n - {\left(C_{2} t + C_{3}\right)}^{2 \, n} C}{{\left(C_{2} t + C_{3}\right)}^{n} {\left(2 \, w + 3\right)}} $$

which is of the general form :

DEg=diff(H(t=t), t)==sum([(K:=var("k", n=3))[u]*H(t=t)^u for u in range(3)])

$$ \frac{\partial}{\partial t}H\left(t\right) = k_{2} H\left(t\right)^{2} + k_{1} H\left(t\right) + k_{0} $$

which can be (awkwardly) solved by Sage :

sage: desolve(DEg, H(t=t), ivar=t)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)

[Snip... ]

TypeError: Computation failed since Maxima requested additional constraints; using the 'assume' command before evaluation *may* help (example of legal syntax is 'assume(4*k0*k2-k1^2>0)', see `assume?` for more details)
Is 4*k0*k2-k1^2 positive or negative?

sage: with assuming (4*k0*k2-k1^2>0): Solp=desolve(DEg, H(t=t), ivar=t) ; Solp
2*arctan((2*k2*H(t) + k1)/sqrt(-k1^2 + 4*k0*k2))/sqrt(-k1^2 + 4*k0*k2) == _C + t

An implicit solution, which can be made explicit :

sage: Solp.solve(H(t=t))
[H(t) == 1/2*(sqrt(-k1^2 + 4*k0*k2)*tan(1/2*sqrt(-k1^2 + 4*k0*k2)*_C + 1/2*sqrt(-k1^2 + 4*k0*k2)*t) - k1)/k2]

Similarly :

sage: with assuming (4*k0*k2-k1^2<0): Soln=desolve(DEg, H(t=t), ivar=t).solve(H(t=t)) ; Soln
[H(t) == -1/2*((k1 + sqrt(k1^2 - 4*k0*k2))*e^(sqrt(k1^2 - 4*k0*k2)*_C + sqrt(k1^2 - 4*k0*k2)*t) - k1 + sqrt(k1^2 - 4*k0*k2))/(k2*e^(sqrt(k1^2 - 4*k0*k2)*_C + sqrt(k1^2 - 4*k0*k2)*t) - k2)]
sage: with assuming (4*k0*k2-k1^2==0): Solz=desolve(DEg, H(t=t), ivar=t).solve(H(t=t)) ; Solz
[H(t) == -1/2*(_C*k1 + k1*t + 2)/(_C*k2 + k2*t)]

There you have a solution for the general form of the equation. All that remains to be done is to replace the coefficients by their values (which I am too lazy to do...) and to frame the three solution cases in a cases expression.

EDIT : Nontwithstanding my laziness, here's the abstract of my solution :

D1=dict(zip((H(t=t), sigma(t=t), P(t=t)), var("f", n=3)))
Sol1=solve([des[u].subs(D1) for u in (0, 2)], (f1, f2))
ID1={D1[u]:u for u in D1.keys()}
DE=des[1].subs(D1).subs(Sol1[0]).subs(ID1)
DE1=((-DE)+diff(H(t=t), t)).lhs().subs(D1)
D2={(K:=var("k", n=3))[u[1]]:u[0].simplify_full().factor()
    for u in DE1.coefficients(f0)}
DE2=diff(H(t=t), t)==sum([u[0]*H(t=t)^u[1] for u in DE1.coefficients(f0)])
DEg=diff(H(t=t), t)==sum([(K:=var("k", n=3))[u]*H(t=t)^u for u in range(3)])
from _operator import lt, eq, gt
Sol=[]
for cond in [f(4*k0*k2-k1^2,0) for f in (gt, eq, lt)]:
    with assuming(cond): Sol+=[(cond, desolve(DEg, H, ivar=t).solve(H)[0].rhs().subs(D2).full_simplify().factor())]
Sol2=H(t=t)==cases(Sol2)

which is way yooo wide to LaTeX. Asking Mathematica a bit of help :

Sol2=[]
for cond in [f(4*k0*k2-k1^2,0) for f in (gt, eq, lt)]:
    with assuming(cond): Sol2+=[[cond, desolve(DEg, H, ivar=t).solve(H)[0].rhs().subs(D2)._mathematica_().FullSimplify().sage(locals={("Tan", 1):tan, ("Coth", 1):coth})]]
# Refining display
Sol2[2][1]=Sol2[2][1]._mathematica_().Factor().sage(locals={("Tan", 1):tan, ("Coth", 1):coth})
Sol2=H==cases(Sol2)

Since this site's Mathjax desn't seem to support cases, here is its argument :

$$ \left[-k_{1}^{2} + 4 \, k_{0} k_{2} > 0, -\frac{1}{6} \, \sqrt{-\frac{C_{2}^{2} n^{2}}{{\left(C_{2} t + C_{3}\right)}^{2}} - \frac{6 \, {\left(2 \, {\left(C_{2} t + C_{3}\right)}^{n} C_{2} n + {\left(C_{2} t + C_{3}\right)} {\left(2 \, w + 1\right)}\right)} {\left(C_{2} t + C_{3}\right)}^{n - 1} C}{2 \, w + 3}} \tan\left(\frac{1}{2} \, \sqrt{-\frac{C_{2}^{2} n^{2}}{{\left(C_{2} t + C_{3}\right)}^{2}} - \frac{6 \, {\left(2 \, {\left(C_{2} t + C_{3}\right)}^{n} C_{2} n + {\left(C_{2} t + C_{3}\right)} {\left(2 \, w + 1\right)}\right)} {\left(C_{2} t + C_{3}\right)}^{n - 1} C}{2 \, w + 3}} {\left(C + t\right)}\right) - \frac{C_{2} n}{6 \, {\left(C_{2} t + C_{3}\right)}}\right] $$

$$ \left[-k_{1}^{2} + 4 \, k_{0} k_{2} = 0, -\frac{C_{2} n}{6 \, {\left(C_{2} t + C_{3}\right)}} + \frac{1}{3 \, {\left(C + t\right)}}\right] $$

$$ \left[-k_{1}^{2} + 4 \, k_{0} k_{2} < 0, -\frac{C_{2} n \sqrt{2 \, w + 3} - \sqrt{C_{2}^{2} n^{2} {\left(2 \, w + 3\right)} + 12 \, {\left(C_{2} t + C_{3}\right)}^{2 \, n + 1} C C_{2} n + 6 \, {\left(C_{2} t + C_{3}\right)}^{n + 2} C {\left(2 \, w + 1\right)}} \coth\left(\frac{\sqrt{C_{2}^{2} n^{2} {\left(2 \, w + 3\right)} + 12 \, {\left(C_{2} t + C_{3}\right)}^{2 \, n + 1} C C_{2} n + 6 \, {\left(C_{2} t + C_{3}\right)}^{n + 2} C {\left(2 \, w + 1\right)}} {\left(C + t\right)}}{2 \, {\left(C_{2} t + C_{3}\right)} \sqrt{2 \, w + 3}}\right)}{6 \, {\left(C_{2} t + C_{3}\right)} \sqrt{2 \, w + 3}}\right] $$

HTH,

edit flag offensive delete link more

Comments

This is a realy great answer, thank you very much, it is appreciated! I have some other related questions. In the begining, all (or most) equations were differential equations, but I tried many ansatz to try to simplify, which gave this mess of constants.

So If I try the general system of equation, where H, \phi, \Sigma and P (P is here only to not have a second order DE, the documentation said desolve was only for first order DE) are all functions of "t", only V(\phi) is fixed, I get solutions in the form of inverse Laplace transform. I get this even with initial conditions (I'm not sure if I'm using them properly) (with desolve_system):

ics = [oo, 0,  0,  0, 0]

Is there a way to fix the boundary conditions to solve these inverse Laplace transforms? Do I have to add more ics?

roux_de_secours gravatar imageroux_de_secours ( 2024-02-21 22:42:35 +0100 )edit

@roux_de_secours : Could we have the original statement of the problem ?I. e. the original system of differential equations (whatever their order) ?

Is there a way to fix the boundary conditions to solve these inverse Laplace transforms?

Probably not : the ilt in the answer means that Maxima is unable to find an explicit expression of this transform. The problem is not with the boundary conditions, but with the solution(s) (it|them)sel(f|ves). But, as we saw, Sage can call orger algorithm leading to a different expression of the solution(s).

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2024-02-23 13:23:14 +0100 )edit

The original equations are

t = var('t')
var('w C1 C2 C3 n ')
H = function('H')(t)
phi = function('\\phi')(t) #(C2 * t + C3)^n
sigma = function('\\Sigma')(t)
V = C1*phi^2 #function('V')(phi)

de1 = H^2-w/6 * (diff(phi,t)/phi)^2 - V/(6*phi) - 2*sigma/3 + H*diff(phi, t)/phi
de2 = diff(H, t) + w/6 * (diff(phi,t)/phi)^2 - V/(3*phi) + 2*sigma/3 + 2*H^2 + (diff(phi, t, 2) + 3*H*diff(phi, t))/(2*phi)
de3 = diff(phi, t, 2) + 3*H*diff(phi, t) + (phi * diff(V, t) - 2*V)/(2*w+3)

Where, originally, \phi and V(\phi) are undefined, but I used the ansatz \phi = (at+b)^n and V = C/2 \phi^2 to have less unknown variables, and for physical reasons. I'm not so sure how many unknowns my equations permit me ...(more)

roux_de_secours gravatar imageroux_de_secours ( 2024-02-23 16:50:03 +0100 )edit

The original system (with phi unknown) cnnot be solved by Sagemath (nor by Sympy nor Mathematica, BTW). giving phi a value (as you did in your initial question) allows to solve. This suggests that there is not enough information in the original system to solve it.

I do not (yet) understand why...

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2024-02-23 23:01:58 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2024-02-20 18:19:47 +0100

Seen: 250 times

Last updated: Feb 21