sympyで「連立微分方程式(Simultaneous differential equation)」をやってみた

何かいい方法があれば教えて下さい。

①やさしい連立微分方程式＞1.1 具体的に解いてみよう

fullscript.py

from sympy import *
var('C1 C2')
t = symbols('t')
x1, x2 = symbols('x1, x2', cls=Function)
eq = (Eq(x1(t).diff(t),4*x1(t)-2*x2(t)), Eq(x2(t).diff(t),x1(t)+x2(t)))
ans=dsolve(eq)
f=ans[0].rhs.subs( C1,-C1/2).subs( C2,-C2)
g=ans[1].rhs.subs( C1,-C1/2).subs( C2,-C2)
print("#x1=",f)
print("#x2=",g)
#x1= C1*exp(2*t) + 2*C2*exp(3*t)
#x2= C1*exp(2*t) + C2*exp(3*t)
# in dsolve raise NotImplementedError
NotImplementedError

②線形代数II/連立線形微分方程式＞1階連立線形微分方程式

fullscript.py

from sympy import *
var('a b c d')
t = symbols('t')
x1, x2 = symbols('x1, x2', cls=Function)
eq = (Eq(x1(t).diff(t),a*x1(t)+b*x2(t)), Eq(x2(t).diff(t),c*x1(t)+d*x2(t)))
ans=dsolve(eq)
print(ans)
# TypeError: cannot determine truth value of Relational

fullscript.py

from sympy import *
var('a b c d')
t = symbols('t')
x1, x2 = symbols('x1, x2', cls=Function)
var('C1 C2')
a=3
b=1
c=1
d=3
eq = (Eq(x1(t).diff(t),a*x1(t)+b*x2(t)), Eq(x2(t).diff(t),c*x1(t)+d*x2(t)))
ans=dsolve(eq)
# print(ans)
f=ans[0].rhs
g=ans[1].rhs
dict=solve([Eq(f.subs(t,0),6), Eq(g.subs(t,0),2)], [C1, C2])
f=f.subs( [ ( C1, dict[C1] ), ( C2, dict[C2] ) ] )
g=g.subs( [ ( C1, dict[C1] ), ( C2, dict[C2] ) ] )
print("#f=",f)
print("#g=",g)
#f= 4*exp(4*t) + 2*exp(2*t)
#g= 4*exp(4*t) - 2*exp(2*t)

③線形代数II/連立線形微分方程式＞2階連立線形微分方程式

fullscript.py

from sympy import *
var('a b c d')
t = symbols('t')
k = symbols('k')
m = symbols('m')
x1, x2 = symbols('x1, x2', cls=Function)
var('C1 C2')
eq = (Eq(m*diff(x1(t),t,2),-k*x1(t)-k*(x1(t)-x2(t))),Eq(m*diff(x2(t),t,2),-k*x2(t)-k*(x2(t)-x1(t))) )
ans=dsolve(eq)
print(ans)
# sympy.polys.polyerrors.PolynomialError: only univariate polynomials are allowed