Ask Your Question
0

How to solve a system of differential equations

asked 2024-06-04 18:37:41 +0100

Cyrille gravatar image

updated 2024-06-04 18:44:47 +0100

I wonder where is my mistake. This system is linear has only two equations. I want the analytical solution but Sagemath returns an error

var('X Y X0 Y0 t a0 b0 A0 a1 b1 A1')
X = function('X')(t)
Y = function('Y')(t)
de1=(diff(X,t) - a0*X - b0*Y + A0 == 0)
de2=(diff(Y,t) - a1*X - b1*Y + A1 == 0)
show(de1)
show(de2)
sol=desolve_system([de1,de2],[X,Y],[0,X0,Y0],t)
show(sol)
edit retag flag offensive close merge delete

2 Answers

Sort by » oldest newest most voted
1

answered 2024-06-04 19:20:57 +0100

achrzesz gravatar image

Add some assumptions, for example

assume(b1^2-2*a0*b1+4*a1*b0>0)
assume(b1^2-2*a0*b1+4*a1*b0+a0^2>0)
edit flag offensive delete link more
0

answered 2024-06-06 15:40:06 +0100

Emmanuel Charpentier gravatar image

updated 2024-06-08 09:46:26 +0100

Alternative to @achresz's answer : change tools. Running :

X0, Y0, t, a0, b0, A0, a1, b1, A1, = SR.var('X0, Y0, t, a0, b0, A0, a1, b1, A1')
X, Y = function("X, Y")
de1 = X(t).diff(t) == a0*X(t)+b0*Y(t)-A0
de2 = Y(t).diff(t) == a1*X(t)+b1*Y(t)-A1
import sympy
%time foo = sympy.dsolve(*map(sympy.sympify, ([de1, de2], [X(t), Y(t)])))
%time tfoo=[u.rewrite("sin") for u in foo]
%time sfoo=[u._sage_() for u in foo]
%time stfoo=[u._sage_() for u in tfoo]

gives :

CPU times: user 1.78 s, sys: 8.09 ms, total: 1.79 s
Wall time: 1.79 s
CPU times: user 24.9 ms, sys: 56 µs, total: 24.9 ms
Wall time: 24.9 ms
CPU times: user 10 ms, sys: 0 ns, total: 10 ms
Wall time: 10 ms
CPU times: user 7.47 ms, sys: 3.89 ms, total: 11.4 ms
Wall time: 11.4 ms

The expressions are somewhat bulky, and neither Sympy nor Sage can lighten them significatively.

In this case, "the competition" (the gratis-but-not-free Wolfram Engine) does better. Running :

M1=mathematica("Sys={X'[t]==a0*X[t]+b0*Y[t]-A0, Y'[t]==a0*X[t]+b1*Y[t]-A1}")
M2=mathematica("TSol=Timing[FullSimplify[DSolve[Sys, {X[t], Y[t]}, t]]]")

gives

CPU times: user 215 µs, sys: 13 µs, total: 228 µs
Wall time: 12.2 s

(half a magnitude slower), but :

sage: M2[2]
{{X[t] -> (A1*b0 - A0*b1)/(a0*b0 - a0*b1) + E^(((a0 + b1)*t)/2)*
     (C[1]*Cosh[(Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]*t)/2] + 
      (((a0 - b1)*C[1] + 2*b0*C[2])*
        Sinh[(Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]*t)/2])/
       Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]), 
  Y[t] -> (A0 - A1)/(b0 - b1) + E^(((a0 + b1)*t)/2)*
     (C[2]*Cosh[(Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]*t)/2] + 
      ((2*a0*C[1] - a0*C[2] + b1*C[2])*
        Sinh[(Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]*t)/2])/
       Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2])}}

waaaay lighter...

Checking the identity of the solutions symbolically is problematic : (M2[2][1][1][2]-stfoo[0].rhs()._mathematica_()).FullSimplify() doesn't return after several minutes.

Checking them numerically is (lazily) left as an exercise to the reader.

EDIT : Fixed a typo (mis-pasting) in Mathematica's code.

HTH,

edit flag offensive delete link more

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-06-04 18:37:41 +0100

Seen: 242 times

Last updated: Jun 08