test_ode.py 47 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105
  1. from sympy.core.function import (Derivative, Function, Subs, diff)
  2. from sympy.core.numbers import (E, I, Rational, pi)
  3. from sympy.core.relational import Eq
  4. from sympy.core.singleton import S
  5. from sympy.core.symbol import (Symbol, symbols)
  6. from sympy.functions.elementary.complexes import (im, re)
  7. from sympy.functions.elementary.exponential import (exp, log)
  8. from sympy.functions.elementary.hyperbolic import acosh
  9. from sympy.functions.elementary.miscellaneous import sqrt
  10. from sympy.functions.elementary.trigonometric import (atan2, cos, sin, tan)
  11. from sympy.integrals.integrals import Integral
  12. from sympy.polys.polytools import Poly
  13. from sympy.series.order import O
  14. from sympy.simplify.radsimp import collect
  15. from sympy.solvers.ode import (classify_ode,
  16. homogeneous_order, dsolve)
  17. from sympy.solvers.ode.subscheck import checkodesol
  18. from sympy.solvers.ode.ode import (classify_sysode,
  19. constant_renumber, constantsimp, get_numbered_constants, solve_ics)
  20. from sympy.solvers.ode.nonhomogeneous import _undetermined_coefficients_match
  21. from sympy.solvers.ode.single import LinearCoefficients
  22. from sympy.solvers.deutils import ode_order
  23. from sympy.testing.pytest import XFAIL, raises, slow, SKIP
  24. from sympy.utilities.misc import filldedent
  25. C0, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10 = symbols('C0:11')
  26. u, x, y, z = symbols('u,x:z', real=True)
  27. f = Function('f')
  28. g = Function('g')
  29. h = Function('h')
  30. # Note: Examples which were specifically testing Single ODE solver are moved to test_single.py
  31. # and all the system of ode examples are moved to test_systems.py
  32. # Note: the tests below may fail (but still be correct) if ODE solver,
  33. # the integral engine, solve(), or even simplify() changes. Also, in
  34. # differently formatted solutions, the arbitrary constants might not be
  35. # equal. Using specific hints in tests can help to avoid this.
  36. # Tests of order higher than 1 should run the solutions through
  37. # constant_renumber because it will normalize it (constant_renumber causes
  38. # dsolve() to return different results on different machines)
  39. def test_get_numbered_constants():
  40. with raises(ValueError):
  41. get_numbered_constants(None)
  42. def test_dsolve_all_hint():
  43. eq = f(x).diff(x)
  44. output = dsolve(eq, hint='all')
  45. # Match the Dummy variables:
  46. sol1 = output['separable_Integral']
  47. _y = sol1.lhs.args[1][0]
  48. sol1 = output['1st_homogeneous_coeff_subs_dep_div_indep_Integral']
  49. _u1 = sol1.rhs.args[1].args[1][0]
  50. expected = {'Bernoulli_Integral': Eq(f(x), C1 + Integral(0, x)),
  51. '1st_homogeneous_coeff_best': Eq(f(x), C1),
  52. 'Bernoulli': Eq(f(x), C1),
  53. 'nth_algebraic': Eq(f(x), C1),
  54. 'nth_linear_euler_eq_homogeneous': Eq(f(x), C1),
  55. 'nth_linear_constant_coeff_homogeneous': Eq(f(x), C1),
  56. 'separable': Eq(f(x), C1),
  57. '1st_homogeneous_coeff_subs_indep_div_dep': Eq(f(x), C1),
  58. 'nth_algebraic_Integral': Eq(f(x), C1),
  59. '1st_linear': Eq(f(x), C1),
  60. '1st_linear_Integral': Eq(f(x), C1 + Integral(0, x)),
  61. '1st_exact': Eq(f(x), C1),
  62. '1st_exact_Integral': Eq(Subs(Integral(0, x) + Integral(1, _y), _y, f(x)), C1),
  63. 'lie_group': Eq(f(x), C1),
  64. '1st_homogeneous_coeff_subs_dep_div_indep': Eq(f(x), C1),
  65. '1st_homogeneous_coeff_subs_dep_div_indep_Integral': Eq(log(x), C1 + Integral(-1/_u1, (_u1, f(x)/x))),
  66. '1st_power_series': Eq(f(x), C1),
  67. 'separable_Integral': Eq(Integral(1, (_y, f(x))), C1 + Integral(0, x)),
  68. '1st_homogeneous_coeff_subs_indep_div_dep_Integral': Eq(f(x), C1),
  69. 'best': Eq(f(x), C1),
  70. 'best_hint': 'nth_algebraic',
  71. 'default': 'nth_algebraic',
  72. 'order': 1}
  73. assert output == expected
  74. assert dsolve(eq, hint='best') == Eq(f(x), C1)
  75. def test_dsolve_ics():
  76. # Maybe this should just use one of the solutions instead of raising...
  77. with raises(NotImplementedError):
  78. dsolve(f(x).diff(x) - sqrt(f(x)), ics={f(1):1})
  79. @slow
  80. def test_dsolve_options():
  81. eq = x*f(x).diff(x) + f(x)
  82. a = dsolve(eq, hint='all')
  83. b = dsolve(eq, hint='all', simplify=False)
  84. c = dsolve(eq, hint='all_Integral')
  85. keys = ['1st_exact', '1st_exact_Integral', '1st_homogeneous_coeff_best',
  86. '1st_homogeneous_coeff_subs_dep_div_indep',
  87. '1st_homogeneous_coeff_subs_dep_div_indep_Integral',
  88. '1st_homogeneous_coeff_subs_indep_div_dep',
  89. '1st_homogeneous_coeff_subs_indep_div_dep_Integral', '1st_linear',
  90. '1st_linear_Integral', 'Bernoulli', 'Bernoulli_Integral',
  91. 'almost_linear', 'almost_linear_Integral', 'best', 'best_hint',
  92. 'default', 'factorable', 'lie_group',
  93. 'nth_linear_euler_eq_homogeneous', 'order',
  94. 'separable', 'separable_Integral']
  95. Integral_keys = ['1st_exact_Integral',
  96. '1st_homogeneous_coeff_subs_dep_div_indep_Integral',
  97. '1st_homogeneous_coeff_subs_indep_div_dep_Integral', '1st_linear_Integral',
  98. 'Bernoulli_Integral', 'almost_linear_Integral', 'best', 'best_hint', 'default',
  99. 'factorable', 'nth_linear_euler_eq_homogeneous',
  100. 'order', 'separable_Integral']
  101. assert sorted(a.keys()) == keys
  102. assert a['order'] == ode_order(eq, f(x))
  103. assert a['best'] == Eq(f(x), C1/x)
  104. assert dsolve(eq, hint='best') == Eq(f(x), C1/x)
  105. assert a['default'] == 'factorable'
  106. assert a['best_hint'] == 'factorable'
  107. assert not a['1st_exact'].has(Integral)
  108. assert not a['separable'].has(Integral)
  109. assert not a['1st_homogeneous_coeff_best'].has(Integral)
  110. assert not a['1st_homogeneous_coeff_subs_dep_div_indep'].has(Integral)
  111. assert not a['1st_homogeneous_coeff_subs_indep_div_dep'].has(Integral)
  112. assert not a['1st_linear'].has(Integral)
  113. assert a['1st_linear_Integral'].has(Integral)
  114. assert a['1st_exact_Integral'].has(Integral)
  115. assert a['1st_homogeneous_coeff_subs_dep_div_indep_Integral'].has(Integral)
  116. assert a['1st_homogeneous_coeff_subs_indep_div_dep_Integral'].has(Integral)
  117. assert a['separable_Integral'].has(Integral)
  118. assert sorted(b.keys()) == keys
  119. assert b['order'] == ode_order(eq, f(x))
  120. assert b['best'] == Eq(f(x), C1/x)
  121. assert dsolve(eq, hint='best', simplify=False) == Eq(f(x), C1/x)
  122. assert b['default'] == 'factorable'
  123. assert b['best_hint'] == 'factorable'
  124. assert a['separable'] != b['separable']
  125. assert a['1st_homogeneous_coeff_subs_dep_div_indep'] != \
  126. b['1st_homogeneous_coeff_subs_dep_div_indep']
  127. assert a['1st_homogeneous_coeff_subs_indep_div_dep'] != \
  128. b['1st_homogeneous_coeff_subs_indep_div_dep']
  129. assert not b['1st_exact'].has(Integral)
  130. assert not b['separable'].has(Integral)
  131. assert not b['1st_homogeneous_coeff_best'].has(Integral)
  132. assert not b['1st_homogeneous_coeff_subs_dep_div_indep'].has(Integral)
  133. assert not b['1st_homogeneous_coeff_subs_indep_div_dep'].has(Integral)
  134. assert not b['1st_linear'].has(Integral)
  135. assert b['1st_linear_Integral'].has(Integral)
  136. assert b['1st_exact_Integral'].has(Integral)
  137. assert b['1st_homogeneous_coeff_subs_dep_div_indep_Integral'].has(Integral)
  138. assert b['1st_homogeneous_coeff_subs_indep_div_dep_Integral'].has(Integral)
  139. assert b['separable_Integral'].has(Integral)
  140. assert sorted(c.keys()) == Integral_keys
  141. raises(ValueError, lambda: dsolve(eq, hint='notarealhint'))
  142. raises(ValueError, lambda: dsolve(eq, hint='Liouville'))
  143. assert dsolve(f(x).diff(x) - 1/f(x)**2, hint='all')['best'] == \
  144. dsolve(f(x).diff(x) - 1/f(x)**2, hint='best')
  145. assert dsolve(f(x) + f(x).diff(x) + sin(x).diff(x) + 1, f(x),
  146. hint="1st_linear_Integral") == \
  147. Eq(f(x), (C1 + Integral((-sin(x).diff(x) - 1)*
  148. exp(Integral(1, x)), x))*exp(-Integral(1, x)))
  149. def test_classify_ode():
  150. assert classify_ode(f(x).diff(x, 2), f(x)) == \
  151. (
  152. 'nth_algebraic',
  153. 'nth_linear_constant_coeff_homogeneous',
  154. 'nth_linear_euler_eq_homogeneous',
  155. 'Liouville',
  156. '2nd_power_series_ordinary',
  157. 'nth_algebraic_Integral',
  158. 'Liouville_Integral',
  159. )
  160. assert classify_ode(f(x), f(x)) == ('nth_algebraic', 'nth_algebraic_Integral')
  161. assert classify_ode(Eq(f(x).diff(x), 0), f(x)) == (
  162. 'nth_algebraic',
  163. 'separable',
  164. '1st_exact',
  165. '1st_linear',
  166. 'Bernoulli',
  167. '1st_homogeneous_coeff_best',
  168. '1st_homogeneous_coeff_subs_indep_div_dep',
  169. '1st_homogeneous_coeff_subs_dep_div_indep',
  170. '1st_power_series', 'lie_group',
  171. 'nth_linear_constant_coeff_homogeneous',
  172. 'nth_linear_euler_eq_homogeneous',
  173. 'nth_algebraic_Integral',
  174. 'separable_Integral',
  175. '1st_exact_Integral',
  176. '1st_linear_Integral',
  177. 'Bernoulli_Integral',
  178. '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
  179. '1st_homogeneous_coeff_subs_dep_div_indep_Integral')
  180. assert classify_ode(f(x).diff(x)**2, f(x)) == ('factorable',
  181. 'nth_algebraic',
  182. 'separable',
  183. '1st_exact',
  184. '1st_linear',
  185. 'Bernoulli',
  186. '1st_homogeneous_coeff_best',
  187. '1st_homogeneous_coeff_subs_indep_div_dep',
  188. '1st_homogeneous_coeff_subs_dep_div_indep',
  189. '1st_power_series',
  190. 'lie_group',
  191. 'nth_linear_euler_eq_homogeneous',
  192. 'nth_algebraic_Integral',
  193. 'separable_Integral',
  194. '1st_exact_Integral',
  195. '1st_linear_Integral',
  196. 'Bernoulli_Integral',
  197. '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
  198. '1st_homogeneous_coeff_subs_dep_div_indep_Integral')
  199. # issue 4749: f(x) should be cleared from highest derivative before classifying
  200. a = classify_ode(Eq(f(x).diff(x) + f(x), x), f(x))
  201. b = classify_ode(f(x).diff(x)*f(x) + f(x)*f(x) - x*f(x), f(x))
  202. c = classify_ode(f(x).diff(x)/f(x) + f(x)/f(x) - x/f(x), f(x))
  203. assert a == ('1st_exact',
  204. '1st_linear',
  205. 'Bernoulli',
  206. 'almost_linear',
  207. '1st_power_series', "lie_group",
  208. 'nth_linear_constant_coeff_undetermined_coefficients',
  209. 'nth_linear_constant_coeff_variation_of_parameters',
  210. '1st_exact_Integral',
  211. '1st_linear_Integral',
  212. 'Bernoulli_Integral',
  213. 'almost_linear_Integral',
  214. 'nth_linear_constant_coeff_variation_of_parameters_Integral')
  215. assert b == ('factorable',
  216. '1st_linear',
  217. 'Bernoulli',
  218. '1st_power_series',
  219. 'lie_group',
  220. 'nth_linear_constant_coeff_undetermined_coefficients',
  221. 'nth_linear_constant_coeff_variation_of_parameters',
  222. '1st_linear_Integral',
  223. 'Bernoulli_Integral',
  224. 'nth_linear_constant_coeff_variation_of_parameters_Integral')
  225. assert c == ('factorable',
  226. '1st_linear',
  227. 'Bernoulli',
  228. '1st_power_series',
  229. 'lie_group',
  230. 'nth_linear_constant_coeff_undetermined_coefficients',
  231. 'nth_linear_constant_coeff_variation_of_parameters',
  232. '1st_linear_Integral',
  233. 'Bernoulli_Integral',
  234. 'nth_linear_constant_coeff_variation_of_parameters_Integral')
  235. assert classify_ode(
  236. 2*x*f(x)*f(x).diff(x) + (1 + x)*f(x)**2 - exp(x), f(x)
  237. ) == ('factorable', '1st_exact', 'Bernoulli', 'almost_linear', 'lie_group',
  238. '1st_exact_Integral', 'Bernoulli_Integral', 'almost_linear_Integral')
  239. assert 'Riccati_special_minus2' in \
  240. classify_ode(2*f(x).diff(x) + f(x)**2 - f(x)/x + 3*x**(-2), f(x))
  241. raises(ValueError, lambda: classify_ode(x + f(x, y).diff(x).diff(
  242. y), f(x, y)))
  243. # issue 5176
  244. k = Symbol('k')
  245. assert classify_ode(f(x).diff(x)/(k*f(x) + k*x*f(x)) + 2*f(x)/(k*f(x) +
  246. k*x*f(x)) + x*f(x).diff(x)/(k*f(x) + k*x*f(x)) + z, f(x)) == \
  247. ('factorable', 'separable', '1st_exact', '1st_linear', 'Bernoulli',
  248. '1st_power_series', 'lie_group', 'separable_Integral', '1st_exact_Integral',
  249. '1st_linear_Integral', 'Bernoulli_Integral')
  250. # preprocessing
  251. ans = ('factorable', 'nth_algebraic', 'separable', '1st_exact', '1st_linear', 'Bernoulli',
  252. '1st_homogeneous_coeff_best',
  253. '1st_homogeneous_coeff_subs_indep_div_dep',
  254. '1st_homogeneous_coeff_subs_dep_div_indep',
  255. '1st_power_series', 'lie_group',
  256. 'nth_linear_constant_coeff_undetermined_coefficients',
  257. 'nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients',
  258. 'nth_linear_constant_coeff_variation_of_parameters',
  259. 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters',
  260. 'nth_algebraic_Integral',
  261. 'separable_Integral', '1st_exact_Integral',
  262. '1st_linear_Integral',
  263. 'Bernoulli_Integral',
  264. '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
  265. '1st_homogeneous_coeff_subs_dep_div_indep_Integral',
  266. 'nth_linear_constant_coeff_variation_of_parameters_Integral',
  267. 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral')
  268. # w/o f(x) given
  269. assert classify_ode(diff(f(x) + x, x) + diff(f(x), x)) == ans
  270. # w/ f(x) and prep=True
  271. assert classify_ode(diff(f(x) + x, x) + diff(f(x), x), f(x),
  272. prep=True) == ans
  273. assert classify_ode(Eq(2*x**3*f(x).diff(x), 0), f(x)) == \
  274. ('factorable', 'nth_algebraic', 'separable', '1st_exact',
  275. '1st_linear', 'Bernoulli', '1st_power_series',
  276. 'lie_group', 'nth_linear_euler_eq_homogeneous',
  277. 'nth_algebraic_Integral', 'separable_Integral', '1st_exact_Integral',
  278. '1st_linear_Integral', 'Bernoulli_Integral')
  279. assert classify_ode(Eq(2*f(x)**3*f(x).diff(x), 0), f(x)) == \
  280. ('factorable', 'nth_algebraic', 'separable', '1st_exact', '1st_linear',
  281. 'Bernoulli', '1st_power_series', 'lie_group', 'nth_algebraic_Integral',
  282. 'separable_Integral', '1st_exact_Integral', '1st_linear_Integral',
  283. 'Bernoulli_Integral')
  284. # test issue 13864
  285. assert classify_ode(Eq(diff(f(x), x) - f(x)**x, 0), f(x)) == \
  286. ('1st_power_series', 'lie_group')
  287. assert isinstance(classify_ode(Eq(f(x), 5), f(x), dict=True), dict)
  288. #This is for new behavior of classify_ode when called internally with default, It should
  289. # return the first hint which matches therefore, 'ordered_hints' key will not be there.
  290. assert sorted(classify_ode(Eq(f(x).diff(x), 0), f(x), dict=True).keys()) == \
  291. ['default', 'nth_linear_constant_coeff_homogeneous', 'order']
  292. a = classify_ode(2*x*f(x)*f(x).diff(x) + (1 + x)*f(x)**2 - exp(x), f(x), dict=True, hint='Bernoulli')
  293. assert sorted(a.keys()) == ['Bernoulli', 'Bernoulli_Integral', 'default', 'order', 'ordered_hints']
  294. # test issue 22155
  295. a = classify_ode(f(x).diff(x) - exp(f(x) - x), f(x))
  296. assert a == ('separable',
  297. '1st_exact', '1st_power_series',
  298. 'lie_group', 'separable_Integral',
  299. '1st_exact_Integral')
  300. def test_classify_ode_ics():
  301. # Dummy
  302. eq = f(x).diff(x, x) - f(x)
  303. # Not f(0) or f'(0)
  304. ics = {x: 1}
  305. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  306. ############################
  307. # f(0) type (AppliedUndef) #
  308. ############################
  309. # Wrong function
  310. ics = {g(0): 1}
  311. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  312. # Contains x
  313. ics = {f(x): 1}
  314. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  315. # Too many args
  316. ics = {f(0, 0): 1}
  317. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  318. # point contains x
  319. ics = {f(0): f(x)}
  320. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  321. # Does not raise
  322. ics = {f(0): f(0)}
  323. classify_ode(eq, f(x), ics=ics)
  324. # Does not raise
  325. ics = {f(0): 1}
  326. classify_ode(eq, f(x), ics=ics)
  327. #####################
  328. # f'(0) type (Subs) #
  329. #####################
  330. # Wrong function
  331. ics = {g(x).diff(x).subs(x, 0): 1}
  332. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  333. # Contains x
  334. ics = {f(y).diff(y).subs(y, x): 1}
  335. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  336. # Wrong variable
  337. ics = {f(y).diff(y).subs(y, 0): 1}
  338. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  339. # Too many args
  340. ics = {f(x, y).diff(x).subs(x, 0): 1}
  341. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  342. # Derivative wrt wrong vars
  343. ics = {Derivative(f(x), x, y).subs(x, 0): 1}
  344. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  345. # point contains x
  346. ics = {f(x).diff(x).subs(x, 0): f(x)}
  347. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  348. # Does not raise
  349. ics = {f(x).diff(x).subs(x, 0): f(x).diff(x).subs(x, 0)}
  350. classify_ode(eq, f(x), ics=ics)
  351. # Does not raise
  352. ics = {f(x).diff(x).subs(x, 0): 1}
  353. classify_ode(eq, f(x), ics=ics)
  354. ###########################
  355. # f'(y) type (Derivative) #
  356. ###########################
  357. # Wrong function
  358. ics = {g(x).diff(x).subs(x, y): 1}
  359. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  360. # Contains x
  361. ics = {f(y).diff(y).subs(y, x): 1}
  362. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  363. # Too many args
  364. ics = {f(x, y).diff(x).subs(x, y): 1}
  365. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  366. # Derivative wrt wrong vars
  367. ics = {Derivative(f(x), x, z).subs(x, y): 1}
  368. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  369. # point contains x
  370. ics = {f(x).diff(x).subs(x, y): f(x)}
  371. raises(ValueError, lambda: classify_ode(eq, f(x), ics=ics))
  372. # Does not raise
  373. ics = {f(x).diff(x).subs(x, 0): f(0)}
  374. classify_ode(eq, f(x), ics=ics)
  375. # Does not raise
  376. ics = {f(x).diff(x).subs(x, y): 1}
  377. classify_ode(eq, f(x), ics=ics)
  378. def test_classify_sysode():
  379. # Here x is assumed to be x(t) and y as y(t) for simplicity.
  380. # Similarly diff(x,t) and diff(y,y) is assumed to be x1 and y1 respectively.
  381. k, l, m, n = symbols('k, l, m, n', Integer=True)
  382. k1, k2, k3, l1, l2, l3, m1, m2, m3 = symbols('k1, k2, k3, l1, l2, l3, m1, m2, m3', Integer=True)
  383. P, Q, R, p, q, r = symbols('P, Q, R, p, q, r', cls=Function)
  384. P1, P2, P3, Q1, Q2, R1, R2 = symbols('P1, P2, P3, Q1, Q2, R1, R2', cls=Function)
  385. x, y, z = symbols('x, y, z', cls=Function)
  386. t = symbols('t')
  387. x1 = diff(x(t),t)
  388. y1 = diff(y(t),t)
  389. eq6 = (Eq(x1, exp(k*x(t))*P(x(t),y(t))), Eq(y1,r(y(t))*P(x(t),y(t))))
  390. sol6 = {'no_of_equation': 2, 'func_coeff': {(0, x(t), 0): 0, (1, x(t), 1): 0, (0, x(t), 1): 1, (1, y(t), 0): 0, \
  391. (1, x(t), 0): 0, (0, y(t), 1): 0, (0, y(t), 0): 0, (1, y(t), 1): 1}, 'type_of_equation': 'type2', 'func': \
  392. [x(t), y(t)], 'is_linear': False, 'eq': [-P(x(t), y(t))*exp(k*x(t)) + Derivative(x(t), t), -P(x(t), \
  393. y(t))*r(y(t)) + Derivative(y(t), t)], 'order': {y(t): 1, x(t): 1}}
  394. assert classify_sysode(eq6) == sol6
  395. eq7 = (Eq(x1, x(t)**2+y(t)/x(t)), Eq(y1, x(t)/y(t)))
  396. sol7 = {'no_of_equation': 2, 'func_coeff': {(0, x(t), 0): 0, (1, x(t), 1): 0, (0, x(t), 1): 1, (1, y(t), 0): 0, \
  397. (1, x(t), 0): -1/y(t), (0, y(t), 1): 0, (0, y(t), 0): -1/x(t), (1, y(t), 1): 1}, 'type_of_equation': 'type3', \
  398. 'func': [x(t), y(t)], 'is_linear': False, 'eq': [-x(t)**2 + Derivative(x(t), t) - y(t)/x(t), -x(t)/y(t) + \
  399. Derivative(y(t), t)], 'order': {y(t): 1, x(t): 1}}
  400. assert classify_sysode(eq7) == sol7
  401. eq8 = (Eq(x1, P1(x(t))*Q1(y(t))*R(x(t),y(t),t)), Eq(y1, P1(x(t))*Q1(y(t))*R(x(t),y(t),t)))
  402. sol8 = {'func': [x(t), y(t)], 'is_linear': False, 'type_of_equation': 'type4', 'eq': \
  403. [-P1(x(t))*Q1(y(t))*R(x(t), y(t), t) + Derivative(x(t), t), -P1(x(t))*Q1(y(t))*R(x(t), y(t), t) + \
  404. Derivative(y(t), t)], 'func_coeff': {(0, y(t), 1): 0, (1, y(t), 1): 1, (1, x(t), 1): 0, (0, y(t), 0): 0, \
  405. (1, x(t), 0): 0, (0, x(t), 0): 0, (1, y(t), 0): 0, (0, x(t), 1): 1}, 'order': {y(t): 1, x(t): 1}, 'no_of_equation': 2}
  406. assert classify_sysode(eq8) == sol8
  407. eq11 = (Eq(x1,x(t)*y(t)**3), Eq(y1,y(t)**5))
  408. sol11 = {'no_of_equation': 2, 'func_coeff': {(0, x(t), 0): -y(t)**3, (1, x(t), 1): 0, (0, x(t), 1): 1, \
  409. (1, y(t), 0): 0, (1, x(t), 0): 0, (0, y(t), 1): 0, (0, y(t), 0): 0, (1, y(t), 1): 1}, 'type_of_equation': \
  410. 'type1', 'func': [x(t), y(t)], 'is_linear': False, 'eq': [-x(t)*y(t)**3 + Derivative(x(t), t), \
  411. -y(t)**5 + Derivative(y(t), t)], 'order': {y(t): 1, x(t): 1}}
  412. assert classify_sysode(eq11) == sol11
  413. eq13 = (Eq(x1,x(t)*y(t)*sin(t)**2), Eq(y1,y(t)**2*sin(t)**2))
  414. sol13 = {'no_of_equation': 2, 'func_coeff': {(0, x(t), 0): -y(t)*sin(t)**2, (1, x(t), 1): 0, (0, x(t), 1): 1, \
  415. (1, y(t), 0): 0, (1, x(t), 0): 0, (0, y(t), 1): 0, (0, y(t), 0): -x(t)*sin(t)**2, (1, y(t), 1): 1}, \
  416. 'type_of_equation': 'type4', 'func': [x(t), y(t)], 'is_linear': False, 'eq': [-x(t)*y(t)*sin(t)**2 + \
  417. Derivative(x(t), t), -y(t)**2*sin(t)**2 + Derivative(y(t), t)], 'order': {y(t): 1, x(t): 1}}
  418. assert classify_sysode(eq13) == sol13
  419. def test_solve_ics():
  420. # Basic tests that things work from dsolve.
  421. assert dsolve(f(x).diff(x) - 1/f(x), f(x), ics={f(1): 2}) == \
  422. Eq(f(x), sqrt(2 * x + 2))
  423. assert dsolve(f(x).diff(x) - f(x), f(x), ics={f(0): 1}) == Eq(f(x), exp(x))
  424. assert dsolve(f(x).diff(x) - f(x), f(x), ics={f(x).diff(x).subs(x, 0): 1}) == Eq(f(x), exp(x))
  425. assert dsolve(f(x).diff(x, x) + f(x), f(x), ics={f(0): 1,
  426. f(x).diff(x).subs(x, 0): 1}) == Eq(f(x), sin(x) + cos(x))
  427. assert dsolve([f(x).diff(x) - f(x) + g(x), g(x).diff(x) - g(x) - f(x)],
  428. [f(x), g(x)], ics={f(0): 1, g(0): 0}) == [Eq(f(x), exp(x)*cos(x)), Eq(g(x), exp(x)*sin(x))]
  429. # Test cases where dsolve returns two solutions.
  430. eq = (x**2*f(x)**2 - x).diff(x)
  431. assert dsolve(eq, f(x), ics={f(1): 0}) == [Eq(f(x),
  432. -sqrt(x - 1)/x), Eq(f(x), sqrt(x - 1)/x)]
  433. assert dsolve(eq, f(x), ics={f(x).diff(x).subs(x, 1): 0}) == [Eq(f(x),
  434. -sqrt(x - S.Half)/x), Eq(f(x), sqrt(x - S.Half)/x)]
  435. eq = cos(f(x)) - (x*sin(f(x)) - f(x)**2)*f(x).diff(x)
  436. assert dsolve(eq, f(x),
  437. ics={f(0):1}, hint='1st_exact', simplify=False) == Eq(x*cos(f(x)) + f(x)**3/3, Rational(1, 3))
  438. assert dsolve(eq, f(x),
  439. ics={f(0):1}, hint='1st_exact', simplify=True) == Eq(x*cos(f(x)) + f(x)**3/3, Rational(1, 3))
  440. assert solve_ics([Eq(f(x), C1*exp(x))], [f(x)], [C1], {f(0): 1}) == {C1: 1}
  441. assert solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2],
  442. {f(0): 1, f(pi/2): 1}) == {C1: 1, C2: 1}
  443. assert solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2],
  444. {f(0): 1, f(x).diff(x).subs(x, 0): 1}) == {C1: 1, C2: 1}
  445. assert solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1}) == \
  446. {C2: 1}
  447. # Some more complicated tests Refer to PR #16098
  448. assert set(dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0, f(x).diff(x).subs(x, 1):0})) == \
  449. {Eq(f(x), 0), Eq(f(x), x ** 3 / 6 - x / 2)}
  450. assert set(dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0})) == \
  451. {Eq(f(x), 0), Eq(f(x), C2*x + x**3/6)}
  452. K, r, f0 = symbols('K r f0')
  453. sol = Eq(f(x), K*f0*exp(r*x)/((-K + f0)*(f0*exp(r*x)/(-K + f0) - 1)))
  454. assert (dsolve(Eq(f(x).diff(x), r * f(x) * (1 - f(x) / K)), f(x), ics={f(0): f0})) == sol
  455. #Order dependent issues Refer to PR #16098
  456. assert set(dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(x).diff(x).subs(x,0):0, f(0):0})) == \
  457. {Eq(f(x), 0), Eq(f(x), x ** 3 / 6)}
  458. assert set(dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0, f(x).diff(x).subs(x,0):0})) == \
  459. {Eq(f(x), 0), Eq(f(x), x ** 3 / 6)}
  460. # XXX: Ought to be ValueError
  461. raises(ValueError, lambda: solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1, f(pi): 1}))
  462. # Degenerate case. f'(0) is identically 0.
  463. raises(ValueError, lambda: solve_ics([Eq(f(x), sqrt(C1 - x**2))], [f(x)], [C1], {f(x).diff(x).subs(x, 0): 0}))
  464. EI, q, L = symbols('EI q L')
  465. # eq = Eq(EI*diff(f(x), x, 4), q)
  466. sols = [Eq(f(x), C1 + C2*x + C3*x**2 + C4*x**3 + q*x**4/(24*EI))]
  467. funcs = [f(x)]
  468. constants = [C1, C2, C3, C4]
  469. # Test both cases, Derivative (the default from f(x).diff(x).subs(x, L)),
  470. # and Subs
  471. ics1 = {f(0): 0,
  472. f(x).diff(x).subs(x, 0): 0,
  473. f(L).diff(L, 2): 0,
  474. f(L).diff(L, 3): 0}
  475. ics2 = {f(0): 0,
  476. f(x).diff(x).subs(x, 0): 0,
  477. Subs(f(x).diff(x, 2), x, L): 0,
  478. Subs(f(x).diff(x, 3), x, L): 0}
  479. solved_constants1 = solve_ics(sols, funcs, constants, ics1)
  480. solved_constants2 = solve_ics(sols, funcs, constants, ics2)
  481. assert solved_constants1 == solved_constants2 == {
  482. C1: 0,
  483. C2: 0,
  484. C3: L**2*q/(4*EI),
  485. C4: -L*q/(6*EI)}
  486. # Allow the ics to refer to f
  487. ics = {f(0): f(0)}
  488. assert dsolve(f(x).diff(x) - f(x), f(x), ics=ics) == Eq(f(x), f(0)*exp(x))
  489. ics = {f(x).diff(x).subs(x, 0): f(x).diff(x).subs(x, 0), f(0): f(0)}
  490. assert dsolve(f(x).diff(x, x) + f(x), f(x), ics=ics) == \
  491. Eq(f(x), f(0)*cos(x) + f(x).diff(x).subs(x, 0)*sin(x))
  492. def test_ode_order():
  493. f = Function('f')
  494. g = Function('g')
  495. x = Symbol('x')
  496. assert ode_order(3*x*exp(f(x)), f(x)) == 0
  497. assert ode_order(x*diff(f(x), x) + 3*x*f(x) - sin(x)/x, f(x)) == 1
  498. assert ode_order(x**2*f(x).diff(x, x) + x*diff(f(x), x) - f(x), f(x)) == 2
  499. assert ode_order(diff(x*exp(f(x)), x, x), f(x)) == 2
  500. assert ode_order(diff(x*diff(x*exp(f(x)), x, x), x), f(x)) == 3
  501. assert ode_order(diff(f(x), x, x), g(x)) == 0
  502. assert ode_order(diff(f(x), x, x)*diff(g(x), x), f(x)) == 2
  503. assert ode_order(diff(f(x), x, x)*diff(g(x), x), g(x)) == 1
  504. assert ode_order(diff(x*diff(x*exp(f(x)), x, x), x), g(x)) == 0
  505. # issue 5835: ode_order has to also work for unevaluated derivatives
  506. # (ie, without using doit()).
  507. assert ode_order(Derivative(x*f(x), x), f(x)) == 1
  508. assert ode_order(x*sin(Derivative(x*f(x)**2, x, x)), f(x)) == 2
  509. assert ode_order(Derivative(x*Derivative(x*exp(f(x)), x, x), x), g(x)) == 0
  510. assert ode_order(Derivative(f(x), x, x), g(x)) == 0
  511. assert ode_order(Derivative(x*exp(f(x)), x, x), f(x)) == 2
  512. assert ode_order(Derivative(f(x), x, x)*Derivative(g(x), x), g(x)) == 1
  513. assert ode_order(Derivative(x*Derivative(f(x), x, x), x), f(x)) == 3
  514. assert ode_order(
  515. x*sin(Derivative(x*Derivative(f(x), x)**2, x, x)), f(x)) == 3
  516. def test_homogeneous_order():
  517. assert homogeneous_order(exp(y/x) + tan(y/x), x, y) == 0
  518. assert homogeneous_order(x**2 + sin(x)*cos(y), x, y) is None
  519. assert homogeneous_order(x - y - x*sin(y/x), x, y) == 1
  520. assert homogeneous_order((x*y + sqrt(x**4 + y**4) + x**2*(log(x) - log(y)))/
  521. (pi*x**Rational(2, 3)*sqrt(y)**3), x, y) == Rational(-1, 6)
  522. assert homogeneous_order(y/x*cos(y/x) - x/y*sin(y/x) + cos(y/x), x, y) == 0
  523. assert homogeneous_order(f(x), x, f(x)) == 1
  524. assert homogeneous_order(f(x)**2, x, f(x)) == 2
  525. assert homogeneous_order(x*y*z, x, y) == 2
  526. assert homogeneous_order(x*y*z, x, y, z) == 3
  527. assert homogeneous_order(x**2*f(x)/sqrt(x**2 + f(x)**2), f(x)) is None
  528. assert homogeneous_order(f(x, y)**2, x, f(x, y), y) == 2
  529. assert homogeneous_order(f(x, y)**2, x, f(x), y) is None
  530. assert homogeneous_order(f(x, y)**2, x, f(x, y)) is None
  531. assert homogeneous_order(f(y, x)**2, x, y, f(x, y)) is None
  532. assert homogeneous_order(f(y), f(x), x) is None
  533. assert homogeneous_order(-f(x)/x + 1/sin(f(x)/ x), f(x), x) == 0
  534. assert homogeneous_order(log(1/y) + log(x**2), x, y) is None
  535. assert homogeneous_order(log(1/y) + log(x), x, y) == 0
  536. assert homogeneous_order(log(x/y), x, y) == 0
  537. assert homogeneous_order(2*log(1/y) + 2*log(x), x, y) == 0
  538. a = Symbol('a')
  539. assert homogeneous_order(a*log(1/y) + a*log(x), x, y) == 0
  540. assert homogeneous_order(f(x).diff(x), x, y) is None
  541. assert homogeneous_order(-f(x).diff(x) + x, x, y) is None
  542. assert homogeneous_order(O(x), x, y) is None
  543. assert homogeneous_order(x + O(x**2), x, y) is None
  544. assert homogeneous_order(x**pi, x) == pi
  545. assert homogeneous_order(x**x, x) is None
  546. raises(ValueError, lambda: homogeneous_order(x*y))
  547. @XFAIL
  548. def test_noncircularized_real_imaginary_parts():
  549. # If this passes, lines numbered 3878-3882 (at the time of this commit)
  550. # of sympy/solvers/ode.py for nth_linear_constant_coeff_homogeneous
  551. # should be removed.
  552. y = sqrt(1+x)
  553. i, r = im(y), re(y)
  554. assert not (i.has(atan2) and r.has(atan2))
  555. def test_collect_respecting_exponentials():
  556. # If this test passes, lines 1306-1311 (at the time of this commit)
  557. # of sympy/solvers/ode.py should be removed.
  558. sol = 1 + exp(x/2)
  559. assert sol == collect( sol, exp(x/3))
  560. def test_undetermined_coefficients_match():
  561. assert _undetermined_coefficients_match(g(x), x) == {'test': False}
  562. assert _undetermined_coefficients_match(sin(2*x + sqrt(5)), x) == \
  563. {'test': True, 'trialset':
  564. {cos(2*x + sqrt(5)), sin(2*x + sqrt(5))}}
  565. assert _undetermined_coefficients_match(sin(x)*cos(x), x) == \
  566. {'test': False}
  567. s = {cos(x), x*cos(x), x**2*cos(x), x**2*sin(x), x*sin(x), sin(x)}
  568. assert _undetermined_coefficients_match(sin(x)*(x**2 + x + 1), x) == \
  569. {'test': True, 'trialset': s}
  570. assert _undetermined_coefficients_match(
  571. sin(x)*x**2 + sin(x)*x + sin(x), x) == {'test': True, 'trialset': s}
  572. assert _undetermined_coefficients_match(
  573. exp(2*x)*sin(x)*(x**2 + x + 1), x
  574. ) == {
  575. 'test': True, 'trialset': {exp(2*x)*sin(x), x**2*exp(2*x)*sin(x),
  576. cos(x)*exp(2*x), x**2*cos(x)*exp(2*x), x*cos(x)*exp(2*x),
  577. x*exp(2*x)*sin(x)}}
  578. assert _undetermined_coefficients_match(1/sin(x), x) == {'test': False}
  579. assert _undetermined_coefficients_match(log(x), x) == {'test': False}
  580. assert _undetermined_coefficients_match(2**(x)*(x**2 + x + 1), x) == \
  581. {'test': True, 'trialset': {2**x, x*2**x, x**2*2**x}}
  582. assert _undetermined_coefficients_match(x**y, x) == {'test': False}
  583. assert _undetermined_coefficients_match(exp(x)*exp(2*x + 1), x) == \
  584. {'test': True, 'trialset': {exp(1 + 3*x)}}
  585. assert _undetermined_coefficients_match(sin(x)*(x**2 + x + 1), x) == \
  586. {'test': True, 'trialset': {x*cos(x), x*sin(x), x**2*cos(x),
  587. x**2*sin(x), cos(x), sin(x)}}
  588. assert _undetermined_coefficients_match(sin(x)*(x + sin(x)), x) == \
  589. {'test': False}
  590. assert _undetermined_coefficients_match(sin(x)*(x + sin(2*x)), x) == \
  591. {'test': False}
  592. assert _undetermined_coefficients_match(sin(x)*tan(x), x) == \
  593. {'test': False}
  594. assert _undetermined_coefficients_match(
  595. x**2*sin(x)*exp(x) + x*sin(x) + x, x
  596. ) == {
  597. 'test': True, 'trialset': {x**2*cos(x)*exp(x), x, cos(x), S.One,
  598. exp(x)*sin(x), sin(x), x*exp(x)*sin(x), x*cos(x), x*cos(x)*exp(x),
  599. x*sin(x), cos(x)*exp(x), x**2*exp(x)*sin(x)}}
  600. assert _undetermined_coefficients_match(4*x*sin(x - 2), x) == {
  601. 'trialset': {x*cos(x - 2), x*sin(x - 2), cos(x - 2), sin(x - 2)},
  602. 'test': True,
  603. }
  604. assert _undetermined_coefficients_match(2**x*x, x) == \
  605. {'test': True, 'trialset': {2**x, x*2**x}}
  606. assert _undetermined_coefficients_match(2**x*exp(2*x), x) == \
  607. {'test': True, 'trialset': {2**x*exp(2*x)}}
  608. assert _undetermined_coefficients_match(exp(-x)/x, x) == \
  609. {'test': False}
  610. # Below are from Ordinary Differential Equations,
  611. # Tenenbaum and Pollard, pg. 231
  612. assert _undetermined_coefficients_match(S(4), x) == \
  613. {'test': True, 'trialset': {S.One}}
  614. assert _undetermined_coefficients_match(12*exp(x), x) == \
  615. {'test': True, 'trialset': {exp(x)}}
  616. assert _undetermined_coefficients_match(exp(I*x), x) == \
  617. {'test': True, 'trialset': {exp(I*x)}}
  618. assert _undetermined_coefficients_match(sin(x), x) == \
  619. {'test': True, 'trialset': {cos(x), sin(x)}}
  620. assert _undetermined_coefficients_match(cos(x), x) == \
  621. {'test': True, 'trialset': {cos(x), sin(x)}}
  622. assert _undetermined_coefficients_match(8 + 6*exp(x) + 2*sin(x), x) == \
  623. {'test': True, 'trialset': {S.One, cos(x), sin(x), exp(x)}}
  624. assert _undetermined_coefficients_match(x**2, x) == \
  625. {'test': True, 'trialset': {S.One, x, x**2}}
  626. assert _undetermined_coefficients_match(9*x*exp(x) + exp(-x), x) == \
  627. {'test': True, 'trialset': {x*exp(x), exp(x), exp(-x)}}
  628. assert _undetermined_coefficients_match(2*exp(2*x)*sin(x), x) == \
  629. {'test': True, 'trialset': {exp(2*x)*sin(x), cos(x)*exp(2*x)}}
  630. assert _undetermined_coefficients_match(x - sin(x), x) == \
  631. {'test': True, 'trialset': {S.One, x, cos(x), sin(x)}}
  632. assert _undetermined_coefficients_match(x**2 + 2*x, x) == \
  633. {'test': True, 'trialset': {S.One, x, x**2}}
  634. assert _undetermined_coefficients_match(4*x*sin(x), x) == \
  635. {'test': True, 'trialset': {x*cos(x), x*sin(x), cos(x), sin(x)}}
  636. assert _undetermined_coefficients_match(x*sin(2*x), x) == \
  637. {'test': True, 'trialset':
  638. {x*cos(2*x), x*sin(2*x), cos(2*x), sin(2*x)}}
  639. assert _undetermined_coefficients_match(x**2*exp(-x), x) == \
  640. {'test': True, 'trialset': {x*exp(-x), x**2*exp(-x), exp(-x)}}
  641. assert _undetermined_coefficients_match(2*exp(-x) - x**2*exp(-x), x) == \
  642. {'test': True, 'trialset': {x*exp(-x), x**2*exp(-x), exp(-x)}}
  643. assert _undetermined_coefficients_match(exp(-2*x) + x**2, x) == \
  644. {'test': True, 'trialset': {S.One, x, x**2, exp(-2*x)}}
  645. assert _undetermined_coefficients_match(x*exp(-x), x) == \
  646. {'test': True, 'trialset': {x*exp(-x), exp(-x)}}
  647. assert _undetermined_coefficients_match(x + exp(2*x), x) == \
  648. {'test': True, 'trialset': {S.One, x, exp(2*x)}}
  649. assert _undetermined_coefficients_match(sin(x) + exp(-x), x) == \
  650. {'test': True, 'trialset': {cos(x), sin(x), exp(-x)}}
  651. assert _undetermined_coefficients_match(exp(x), x) == \
  652. {'test': True, 'trialset': {exp(x)}}
  653. # converted from sin(x)**2
  654. assert _undetermined_coefficients_match(S.Half - cos(2*x)/2, x) == \
  655. {'test': True, 'trialset': {S.One, cos(2*x), sin(2*x)}}
  656. # converted from exp(2*x)*sin(x)**2
  657. assert _undetermined_coefficients_match(
  658. exp(2*x)*(S.Half + cos(2*x)/2), x
  659. ) == {
  660. 'test': True, 'trialset': {exp(2*x)*sin(2*x), cos(2*x)*exp(2*x),
  661. exp(2*x)}}
  662. assert _undetermined_coefficients_match(2*x + sin(x) + cos(x), x) == \
  663. {'test': True, 'trialset': {S.One, x, cos(x), sin(x)}}
  664. # converted from sin(2*x)*sin(x)
  665. assert _undetermined_coefficients_match(cos(x)/2 - cos(3*x)/2, x) == \
  666. {'test': True, 'trialset': {cos(x), cos(3*x), sin(x), sin(3*x)}}
  667. assert _undetermined_coefficients_match(cos(x**2), x) == {'test': False}
  668. assert _undetermined_coefficients_match(2**(x**2), x) == {'test': False}
  669. def test_issue_4785_22462():
  670. from sympy.abc import A
  671. eq = x + A*(x + diff(f(x), x) + f(x)) + diff(f(x), x) + f(x) + 2
  672. assert classify_ode(eq, f(x)) == ('factorable', '1st_exact', '1st_linear',
  673. 'Bernoulli', 'almost_linear', '1st_power_series', 'lie_group',
  674. 'nth_linear_constant_coeff_undetermined_coefficients',
  675. 'nth_linear_constant_coeff_variation_of_parameters',
  676. '1st_exact_Integral', '1st_linear_Integral', 'Bernoulli_Integral',
  677. 'almost_linear_Integral',
  678. 'nth_linear_constant_coeff_variation_of_parameters_Integral')
  679. # issue 4864
  680. eq = (x**2 + f(x)**2)*f(x).diff(x) - 2*x*f(x)
  681. assert classify_ode(eq, f(x)) == ('factorable', '1st_exact',
  682. '1st_homogeneous_coeff_best',
  683. '1st_homogeneous_coeff_subs_indep_div_dep',
  684. '1st_homogeneous_coeff_subs_dep_div_indep',
  685. '1st_power_series',
  686. 'lie_group', '1st_exact_Integral',
  687. '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
  688. '1st_homogeneous_coeff_subs_dep_div_indep_Integral')
  689. def test_issue_4825():
  690. raises(ValueError, lambda: dsolve(f(x, y).diff(x) - y*f(x, y), f(x)))
  691. assert classify_ode(f(x, y).diff(x) - y*f(x, y), f(x), dict=True) == \
  692. {'order': 0, 'default': None, 'ordered_hints': ()}
  693. # See also issue 3793, test Z13.
  694. raises(ValueError, lambda: dsolve(f(x).diff(x), f(y)))
  695. assert classify_ode(f(x).diff(x), f(y), dict=True) == \
  696. {'order': 0, 'default': None, 'ordered_hints': ()}
  697. def test_constant_renumber_order_issue_5308():
  698. from sympy.utilities.iterables import variations
  699. assert constant_renumber(C1*x + C2*y) == \
  700. constant_renumber(C1*y + C2*x) == \
  701. C1*x + C2*y
  702. e = C1*(C2 + x)*(C3 + y)
  703. for a, b, c in variations([C1, C2, C3], 3):
  704. assert constant_renumber(a*(b + x)*(c + y)) == e
  705. def test_constant_renumber():
  706. e1, e2, x, y = symbols("e1:3 x y")
  707. exprs = [e2*x, e1*x + e2*y]
  708. assert constant_renumber(exprs[0]) == e2*x
  709. assert constant_renumber(exprs[0], variables=[x]) == C1*x
  710. assert constant_renumber(exprs[0], variables=[x], newconstants=[C2]) == C2*x
  711. assert constant_renumber(exprs, variables=[x, y]) == [C1*x, C1*y + C2*x]
  712. assert constant_renumber(exprs, variables=[x, y], newconstants=symbols("C3:5")) == [C3*x, C3*y + C4*x]
  713. def test_issue_5770():
  714. k = Symbol("k", real=True)
  715. t = Symbol('t')
  716. w = Function('w')
  717. sol = dsolve(w(t).diff(t, 6) - k**6*w(t), w(t))
  718. assert len([s for s in sol.free_symbols if s.name.startswith('C')]) == 6
  719. assert constantsimp((C1*cos(x) + C2*cos(x))*exp(x), {C1, C2}) == \
  720. C1*cos(x)*exp(x)
  721. assert constantsimp(C1*cos(x) + C2*cos(x) + C3*sin(x), {C1, C2, C3}) == \
  722. C1*cos(x) + C3*sin(x)
  723. assert constantsimp(exp(C1 + x), {C1}) == C1*exp(x)
  724. assert constantsimp(x + C1 + y, {C1, y}) == C1 + x
  725. assert constantsimp(x + C1 + Integral(x, (x, 1, 2)), {C1}) == C1 + x
  726. def test_issue_5112_5430():
  727. assert homogeneous_order(-log(x) + acosh(x), x) is None
  728. assert homogeneous_order(y - log(x), x, y) is None
  729. def test_issue_5095():
  730. f = Function('f')
  731. raises(ValueError, lambda: dsolve(f(x).diff(x)**2, f(x), 'fdsjf'))
  732. def test_homogeneous_function():
  733. f = Function('f')
  734. eq1 = tan(x + f(x))
  735. eq2 = sin((3*x)/(4*f(x)))
  736. eq3 = cos(x*f(x)*Rational(3, 4))
  737. eq4 = log((3*x + 4*f(x))/(5*f(x) + 7*x))
  738. eq5 = exp((2*x**2)/(3*f(x)**2))
  739. eq6 = log((3*x + 4*f(x))/(5*f(x) + 7*x) + exp((2*x**2)/(3*f(x)**2)))
  740. eq7 = sin((3*x)/(5*f(x) + x**2))
  741. assert homogeneous_order(eq1, x, f(x)) == None
  742. assert homogeneous_order(eq2, x, f(x)) == 0
  743. assert homogeneous_order(eq3, x, f(x)) == None
  744. assert homogeneous_order(eq4, x, f(x)) == 0
  745. assert homogeneous_order(eq5, x, f(x)) == 0
  746. assert homogeneous_order(eq6, x, f(x)) == 0
  747. assert homogeneous_order(eq7, x, f(x)) == None
  748. def test_linear_coeff_match():
  749. n, d = z*(2*x + 3*f(x) + 5), z*(7*x + 9*f(x) + 11)
  750. rat = n/d
  751. eq1 = sin(rat) + cos(rat.expand())
  752. obj1 = LinearCoefficients(eq1)
  753. eq2 = rat
  754. obj2 = LinearCoefficients(eq2)
  755. eq3 = log(sin(rat))
  756. obj3 = LinearCoefficients(eq3)
  757. ans = (4, Rational(-13, 3))
  758. assert obj1._linear_coeff_match(eq1, f(x)) == ans
  759. assert obj2._linear_coeff_match(eq2, f(x)) == ans
  760. assert obj3._linear_coeff_match(eq3, f(x)) == ans
  761. # no c
  762. eq4 = (3*x)/f(x)
  763. obj4 = LinearCoefficients(eq4)
  764. # not x and f(x)
  765. eq5 = (3*x + 2)/x
  766. obj5 = LinearCoefficients(eq5)
  767. # denom will be zero
  768. eq6 = (3*x + 2*f(x) + 1)/(3*x + 2*f(x) + 5)
  769. obj6 = LinearCoefficients(eq6)
  770. # not rational coefficient
  771. eq7 = (3*x + 2*f(x) + sqrt(2))/(3*x + 2*f(x) + 5)
  772. obj7 = LinearCoefficients(eq7)
  773. assert obj4._linear_coeff_match(eq4, f(x)) is None
  774. assert obj5._linear_coeff_match(eq5, f(x)) is None
  775. assert obj6._linear_coeff_match(eq6, f(x)) is None
  776. assert obj7._linear_coeff_match(eq7, f(x)) is None
  777. def test_constantsimp_take_problem():
  778. c = exp(C1) + 2
  779. assert len(Poly(constantsimp(exp(C1) + c + c*x, [C1])).gens) == 2
  780. def test_series():
  781. C1 = Symbol("C1")
  782. eq = f(x).diff(x) - f(x)
  783. sol = Eq(f(x), C1 + C1*x + C1*x**2/2 + C1*x**3/6 + C1*x**4/24 +
  784. C1*x**5/120 + O(x**6))
  785. assert dsolve(eq, hint='1st_power_series') == sol
  786. assert checkodesol(eq, sol, order=1)[0]
  787. eq = f(x).diff(x) - x*f(x)
  788. sol = Eq(f(x), C1*x**4/8 + C1*x**2/2 + C1 + O(x**6))
  789. assert dsolve(eq, hint='1st_power_series') == sol
  790. assert checkodesol(eq, sol, order=1)[0]
  791. eq = f(x).diff(x) - sin(x*f(x))
  792. sol = Eq(f(x), (x - 2)**2*(1+ sin(4))*cos(4) + (x - 2)*sin(4) + 2 + O(x**3))
  793. assert dsolve(eq, hint='1st_power_series', ics={f(2): 2}, n=3) == sol
  794. # FIXME: The solution here should be O((x-2)**3) so is incorrect
  795. #assert checkodesol(eq, sol, order=1)[0]
  796. @slow
  797. def test_2nd_power_series_ordinary():
  798. C1, C2 = symbols("C1 C2")
  799. eq = f(x).diff(x, 2) - x*f(x)
  800. assert classify_ode(eq) == ('2nd_linear_airy', '2nd_power_series_ordinary')
  801. sol = Eq(f(x), C2*(x**3/6 + 1) + C1*x*(x**3/12 + 1) + O(x**6))
  802. assert dsolve(eq, hint='2nd_power_series_ordinary') == sol
  803. assert checkodesol(eq, sol) == (True, 0)
  804. sol = Eq(f(x), C2*((x + 2)**4/6 + (x + 2)**3/6 - (x + 2)**2 + 1)
  805. + C1*(x + (x + 2)**4/12 - (x + 2)**3/3 + S(2))
  806. + O(x**6))
  807. assert dsolve(eq, hint='2nd_power_series_ordinary', x0=-2) == sol
  808. # FIXME: Solution should be O((x+2)**6)
  809. # assert checkodesol(eq, sol) == (True, 0)
  810. sol = Eq(f(x), C2*x + C1 + O(x**2))
  811. assert dsolve(eq, hint='2nd_power_series_ordinary', n=2) == sol
  812. assert checkodesol(eq, sol) == (True, 0)
  813. eq = (1 + x**2)*(f(x).diff(x, 2)) + 2*x*(f(x).diff(x)) -2*f(x)
  814. assert classify_ode(eq) == ('factorable', '2nd_hypergeometric', '2nd_hypergeometric_Integral',
  815. '2nd_power_series_ordinary')
  816. sol = Eq(f(x), C2*(-x**4/3 + x**2 + 1) + C1*x + O(x**6))
  817. assert dsolve(eq, hint='2nd_power_series_ordinary') == sol
  818. assert checkodesol(eq, sol) == (True, 0)
  819. eq = f(x).diff(x, 2) + x*(f(x).diff(x)) + f(x)
  820. assert classify_ode(eq) == ('factorable', '2nd_power_series_ordinary',)
  821. sol = Eq(f(x), C2*(x**4/8 - x**2/2 + 1) + C1*x*(-x**2/3 + 1) + O(x**6))
  822. assert dsolve(eq) == sol
  823. # FIXME: checkodesol fails for this solution...
  824. # assert checkodesol(eq, sol) == (True, 0)
  825. eq = f(x).diff(x, 2) + f(x).diff(x) - x*f(x)
  826. assert classify_ode(eq) == ('2nd_power_series_ordinary',)
  827. sol = Eq(f(x), C2*(-x**4/24 + x**3/6 + 1)
  828. + C1*x*(x**3/24 + x**2/6 - x/2 + 1) + O(x**6))
  829. assert dsolve(eq) == sol
  830. # FIXME: checkodesol fails for this solution...
  831. # assert checkodesol(eq, sol) == (True, 0)
  832. eq = f(x).diff(x, 2) + x*f(x)
  833. assert classify_ode(eq) == ('2nd_linear_airy', '2nd_power_series_ordinary')
  834. sol = Eq(f(x), C2*(x**6/180 - x**3/6 + 1) + C1*x*(-x**3/12 + 1) + O(x**7))
  835. assert dsolve(eq, hint='2nd_power_series_ordinary', n=7) == sol
  836. assert checkodesol(eq, sol) == (True, 0)
  837. def test_2nd_power_series_regular():
  838. C1, C2, a = symbols("C1 C2 a")
  839. eq = x**2*(f(x).diff(x, 2)) - 3*x*(f(x).diff(x)) + (4*x + 4)*f(x)
  840. sol = Eq(f(x), C1*x**2*(-16*x**3/9 + 4*x**2 - 4*x + 1) + O(x**6))
  841. assert dsolve(eq, hint='2nd_power_series_regular') == sol
  842. assert checkodesol(eq, sol) == (True, 0)
  843. eq = 4*x**2*(f(x).diff(x, 2)) -8*x**2*(f(x).diff(x)) + (4*x**2 +
  844. 1)*f(x)
  845. sol = Eq(f(x), C1*sqrt(x)*(x**4/24 + x**3/6 + x**2/2 + x + 1) + O(x**6))
  846. assert dsolve(eq, hint='2nd_power_series_regular') == sol
  847. assert checkodesol(eq, sol) == (True, 0)
  848. eq = x**2*(f(x).diff(x, 2)) - x**2*(f(x).diff(x)) + (
  849. x**2 - 2)*f(x)
  850. sol = Eq(f(x), C1*(-x**6/720 - 3*x**5/80 - x**4/8 + x**2/2 + x/2 + 1)/x +
  851. C2*x**2*(-x**3/60 + x**2/20 + x/2 + 1) + O(x**6))
  852. assert dsolve(eq) == sol
  853. assert checkodesol(eq, sol) == (True, 0)
  854. eq = x**2*(f(x).diff(x, 2)) + x*(f(x).diff(x)) + (x**2 - Rational(1, 4))*f(x)
  855. sol = Eq(f(x), C1*(x**4/24 - x**2/2 + 1)/sqrt(x) +
  856. C2*sqrt(x)*(x**4/120 - x**2/6 + 1) + O(x**6))
  857. assert dsolve(eq, hint='2nd_power_series_regular') == sol
  858. assert checkodesol(eq, sol) == (True, 0)
  859. eq = x*f(x).diff(x, 2) + f(x).diff(x) - a*x*f(x)
  860. sol = Eq(f(x), C1*(a**2*x**4/64 + a*x**2/4 + 1) + O(x**6))
  861. assert dsolve(eq, f(x), hint="2nd_power_series_regular") == sol
  862. assert checkodesol(eq, sol) == (True, 0)
  863. eq = f(x).diff(x, 2) + ((1 - x)/x)*f(x).diff(x) + (a/x)*f(x)
  864. sol = Eq(f(x), C1*(-a*x**5*(a - 4)*(a - 3)*(a - 2)*(a - 1)/14400 + \
  865. a*x**4*(a - 3)*(a - 2)*(a - 1)/576 - a*x**3*(a - 2)*(a - 1)/36 + \
  866. a*x**2*(a - 1)/4 - a*x + 1) + O(x**6))
  867. assert dsolve(eq, f(x), hint="2nd_power_series_regular") == sol
  868. assert checkodesol(eq, sol) == (True, 0)
  869. def test_issue_15056():
  870. t = Symbol('t')
  871. C3 = Symbol('C3')
  872. assert get_numbered_constants(Symbol('C1') * Function('C2')(t)) == C3
  873. def test_issue_15913():
  874. eq = -C1/x - 2*x*f(x) - f(x) + Derivative(f(x), x)
  875. sol = C2*exp(x**2 + x) + exp(x**2 + x)*Integral(C1*exp(-x**2 - x)/x, x)
  876. assert checkodesol(eq, sol) == (True, 0)
  877. sol = C1 + C2*exp(-x*y)
  878. eq = Derivative(y*f(x), x) + f(x).diff(x, 2)
  879. assert checkodesol(eq, sol, f(x)) == (True, 0)
  880. def test_issue_16146():
  881. raises(ValueError, lambda: dsolve([f(x).diff(x), g(x).diff(x)], [f(x), g(x), h(x)]))
  882. raises(ValueError, lambda: dsolve([f(x).diff(x), g(x).diff(x)], [f(x)]))
  883. def test_dsolve_remove_redundant_solutions():
  884. eq = (f(x)-2)*f(x).diff(x)
  885. sol = Eq(f(x), C1)
  886. assert dsolve(eq) == sol
  887. eq = (f(x)-sin(x))*(f(x).diff(x, 2))
  888. sol = {Eq(f(x), C1 + C2*x), Eq(f(x), sin(x))}
  889. assert set(dsolve(eq)) == sol
  890. eq = (f(x)**2-2*f(x)+1)*f(x).diff(x, 3)
  891. sol = Eq(f(x), C1 + C2*x + C3*x**2)
  892. assert dsolve(eq) == sol
  893. def test_issue_13060():
  894. A, B = symbols("A B", cls=Function)
  895. t = Symbol("t")
  896. eq = [Eq(Derivative(A(t), t), A(t)*B(t)), Eq(Derivative(B(t), t), A(t)*B(t))]
  897. sol = dsolve(eq)
  898. assert checkodesol(eq, sol) == (True, [0, 0])
  899. def test_issue_22523():
  900. N, s = symbols('N s')
  901. rho = Function('rho')
  902. # intentionally use 4.0 to confirm issue with nfloat
  903. # works here
  904. eqn = 4.0*N*sqrt(N - 1)*rho(s) + (4*s**2*(N - 1) + (N - 2*s*(N - 1))**2
  905. )*Derivative(rho(s), (s, 2))
  906. match = classify_ode(eqn, dict=True, hint='all')
  907. assert match['2nd_power_series_ordinary']['terms'] == 5
  908. C1, C2 = symbols('C1,C2')
  909. sol = dsolve(eqn, hint='2nd_power_series_ordinary')
  910. # there is no r(2.0) in this result
  911. assert filldedent(sol) == filldedent(str('''
  912. Eq(rho(s), C2*(1 - 4.0*s**4*sqrt(N - 1.0)/N + 0.666666666666667*s**4/N
  913. - 2.66666666666667*s**3*sqrt(N - 1.0)/N - 2.0*s**2*sqrt(N - 1.0)/N +
  914. 9.33333333333333*s**4*sqrt(N - 1.0)/N**2 - 0.666666666666667*s**4/N**2
  915. + 2.66666666666667*s**3*sqrt(N - 1.0)/N**2 -
  916. 5.33333333333333*s**4*sqrt(N - 1.0)/N**3) + C1*s*(1.0 -
  917. 1.33333333333333*s**3*sqrt(N - 1.0)/N - 0.666666666666667*s**2*sqrt(N
  918. - 1.0)/N + 1.33333333333333*s**3*sqrt(N - 1.0)/N**2) + O(s**6))'''))
  919. def test_issue_22604():
  920. x1, x2 = symbols('x1, x2', cls = Function)
  921. t, k1, k2, m1, m2 = symbols('t k1 k2 m1 m2', real = True)
  922. k1, k2, m1, m2 = 1, 1, 1, 1
  923. eq1 = Eq(m1*diff(x1(t), t, 2) + k1*x1(t) - k2*(x2(t) - x1(t)), 0)
  924. eq2 = Eq(m2*diff(x2(t), t, 2) + k2*(x2(t) - x1(t)), 0)
  925. eqs = [eq1, eq2]
  926. [x1sol, x2sol] = dsolve(eqs, [x1(t), x2(t)], ics = {x1(0):0, x1(t).diff().subs(t,0):0, \
  927. x2(0):1, x2(t).diff().subs(t,0):0})
  928. assert x1sol == Eq(x1(t), sqrt(3 - sqrt(5))*(sqrt(10) + 5*sqrt(2))*cos(sqrt(2)*t*sqrt(3 - sqrt(5))/2)/20 + \
  929. (-5*sqrt(2) + sqrt(10))*sqrt(sqrt(5) + 3)*cos(sqrt(2)*t*sqrt(sqrt(5) + 3)/2)/20)
  930. assert x2sol == Eq(x2(t), (sqrt(5) + 5)*cos(sqrt(2)*t*sqrt(3 - sqrt(5))/2)/10 + (5 - sqrt(5))*cos(sqrt(2)*t*sqrt(sqrt(5) + 3)/2)/10)
  931. def test_issue_22462():
  932. for de in [
  933. Eq(f(x).diff(x), -20*f(x)**2 - 500*f(x)/7200),
  934. Eq(f(x).diff(x), -2*f(x)**2 - 5*f(x)/7)]:
  935. assert 'Bernoulli' in classify_ode(de, f(x))
  936. def test_issue_23425():
  937. x = symbols('x')
  938. y = Function('y')
  939. eq = Eq(-E**x*y(x).diff().diff() + y(x).diff(), 0)
  940. assert classify_ode(eq) == \
  941. ('Liouville', 'nth_order_reducible', \
  942. '2nd_power_series_ordinary', 'Liouville_Integral')
  943. @SKIP("too slow for @slow")
  944. def test_issue_25820():
  945. x = Symbol('x')
  946. y = Function('y')
  947. eq = y(x)**3*y(x).diff(x, 2) + 49
  948. assert dsolve(eq, y(x)) is not None # doesn't raise