-
-
Notifications
You must be signed in to change notification settings - Fork 17
/
ode.py
7327 lines (6109 loc) · 283 KB
/
ode.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
r"""
This module contains :py:meth:`~diofant.solvers.ode.dsolve` and different helper
functions that it uses.
:py:meth:`~diofant.solvers.ode.dsolve` solves ordinary differential equations.
See the docstring on the various functions for their uses. Note that partial
differential equations support is in ``pde.py``. Note that hint functions
have docstrings describing their various methods, but they are intended for
internal use. Use ``dsolve(ode, func, hint=hint)`` to solve an ODE using a
specific hint. See also the docstring on
:py:meth:`~diofant.solvers.ode.dsolve`.
**Functions in this module**
These are the user functions in this module:
- :py:meth:`~diofant.solvers.ode.dsolve` - Solves ODEs.
- :py:meth:`~diofant.solvers.ode.classify_ode` - Classifies ODEs into
possible hints for :py:meth:`~diofant.solvers.ode.dsolve`.
- :py:meth:`~diofant.solvers.ode.checkodesol` - Checks if an equation is the
solution to an ODE.
- :py:meth:`~diofant.solvers.ode.homogeneous_order` - Returns the
homogeneous order of an expression.
- :py:meth:`~diofant.solvers.ode.infinitesimals` - Returns the infinitesimals
of the Lie group of point transformations of an ODE, such that it is
invariant.
These are the non-solver helper functions that are for internal use. The
user should use the various options to
:py:meth:`~diofant.solvers.ode.dsolve` to obtain the functionality provided
by these functions:
- :py:meth:`~diofant.solvers.ode.odesimp` - Does all forms of ODE
simplification.
- :py:meth:`~diofant.solvers.ode.ode_sol_simplicity` - A key function for
comparing solutions by simplicity.
- :py:meth:`~diofant.solvers.ode.constantsimp` - Simplifies arbitrary
constants.
- :py:meth:`~diofant.solvers.ode.constant_renumber` - Renumber arbitrary
constants.
- :py:meth:`~diofant.solvers.ode._handle_Integral` - Evaluate unevaluated
Integrals.
See also the docstrings of these functions.
**Currently implemented solver methods**
The following methods are implemented for solving ordinary differential
equations. See the docstrings of the various hint functions for more
information on each (run ``help(ode)``):
- 1st order separable differential equations.
- 1st order differential equations whose coefficients or `dx` and `dy` are
functions homogeneous of the same order.
- 1st order exact differential equations.
- 1st order linear differential equations.
- 1st order Bernoulli differential equations.
- Power series solutions for first order differential equations.
- Lie Group method of solving first order differential equations.
- 2nd order Liouville differential equations.
- Power series solutions for second order differential equations
at ordinary and regular singular points.
- `n`\th order linear homogeneous differential equation with constant
coefficients.
- `n`\th order linear inhomogeneous differential equation with constant
coefficients using the method of undetermined coefficients.
- `n`\th order linear inhomogeneous differential equation with constant
coefficients using the method of variation of parameters.
**Philosophy behind this module**
This module is designed to make it easy to add new ODE solving methods without
having to mess with the solving code for other methods. The idea is that
there is a :py:meth:`~diofant.solvers.ode.classify_ode` function, which takes in
an ODE and tells you what hints, if any, will solve the ODE. It does this
without attempting to solve the ODE, so it is fast. Each solving method is a
hint, and it has its own function, named ``ode_<hint>``. That function takes
in the ODE and any match expression gathered by
:py:meth:`~diofant.solvers.ode.classify_ode` and returns a solved result. If
this result has any integrals in it, the hint function will return an
unevaluated :py:class:`~diofant.integrals.integrals.Integral` class.
:py:meth:`~diofant.solvers.ode.dsolve`, which is the user wrapper function
around all of this, will then call :py:meth:`~diofant.solvers.ode.odesimp` on
the result, which, among other things, will attempt to solve the equation for
the dependent variable (the function we are solving for), simplify the
arbitrary constants in the expression, and evaluate any integrals, if the hint
allows it.
**How to add new solution methods**
If you have an ODE that you want :py:meth:`~diofant.solvers.ode.dsolve` to be
able to solve, try to avoid adding special case code here. Instead, try
finding a general method that will solve your ODE, as well as others. This
way, the :py:mod:`~diofant.solvers.ode` module will become more robust, and
unhindered by special case hacks. WolphramAlpha and Maple's
DETools[odeadvisor] function are two resources you can use to classify a
specific ODE. It is also better for a method to work with an `n`\th order ODE
instead of only with specific orders, if possible.
To add a new method, there are a few things that you need to do. First, you
need a hint name for your method. Try to name your hint so that it is
unambiguous with all other methods, including ones that may not be implemented
yet. If your method uses integrals, also include a ``hint_Integral`` hint.
If there is more than one way to solve ODEs with your method, include a hint
for each one, as well as a ``<hint>_best`` hint. Your ``ode_<hint>_best()``
function should choose the best using min with ``ode_sol_simplicity`` as the
key argument. See
:py:meth:`~diofant.solvers.ode.ode_1st_homogeneous_coeff_best`, for example.
The function that uses your method will be called ``ode_<hint>()``, so the
hint must only use characters that are allowed in a Python function name
(alphanumeric characters and the underscore '``_``' character). Include a
function for every hint, except for ``_Integral`` hints
(:py:meth:`~diofant.solvers.ode.dsolve` takes care of those automatically).
Hint names should be all lowercase, unless a word is commonly capitalized
(such as Integral or Bernoulli). If you have a hint that you do not want to
run with ``all_Integral`` that doesn't have an ``_Integral`` counterpart (such
as a best hint that would defeat the purpose of ``all_Integral``), you will
need to remove it manually in the :py:meth:`~diofant.solvers.ode.dsolve` code.
See also the :py:meth:`~diofant.solvers.ode.classify_ode` docstring for
guidelines on writing a hint name.
Determine *in general* how the solutions returned by your method compare with
other methods that can potentially solve the same ODEs. Then, put your hints
in the :py:data:`~diofant.solvers.ode.allhints` tuple in the order that they
should be called. The ordering of this tuple determines which hints are
default. Note that exceptions are ok, because it is easy for the user to
choose individual hints with :py:meth:`~diofant.solvers.ode.dsolve`. In
general, ``_Integral`` variants should go at the end of the list, and
``_best`` variants should go before the various hints they apply to. For
example, the ``undetermined_coefficients`` hint comes before the
``variation_of_parameters`` hint because, even though variation of parameters
is more general than undetermined coefficients, undetermined coefficients
generally returns cleaner results for the ODEs that it can solve than
variation of parameters does, and it does not require integration, so it is
much faster.
Next, you need to have a match expression or a function that matches the type
of the ODE, which you should put in :py:meth:`~diofant.solvers.ode.classify_ode`
(if the match function is more than just a few lines, like
:py:meth:`~diofant.solvers.ode._undetermined_coefficients_match`, it should go
outside of :py:meth:`~diofant.solvers.ode.classify_ode`). It should match the
ODE without solving for it as much as possible, so that
:py:meth:`~diofant.solvers.ode.classify_ode` remains fast and is not hindered by
bugs in solving code. Be sure to consider corner cases. For example, if your
solution method involves dividing by something, make sure you exclude the case
where that division will be 0.
In most cases, the matching of the ODE will also give you the various parts
that you need to solve it. You should put that in a dictionary (``.match()``
will do this for you), and add that as ``matching_hints['hint'] = matchdict``
in the relevant part of :py:meth:`~diofant.solvers.ode.classify_ode`.
:py:meth:`~diofant.solvers.ode.classify_ode` will then send this to
:py:meth:`~diofant.solvers.ode.dsolve`, which will send it to your function as
the ``match`` argument. Your function should be named ``ode_<hint>(eq, func,
order, match)`. If you need to send more information, put it in the ``match``
dictionary. For example, if you had to substitute in a dummy variable in
:py:meth:`~diofant.solvers.ode.classify_ode` to match the ODE, you will need to
pass it to your function using the `match` dict to access it. You can access
the independent variable using ``func.args[0]``, and the dependent variable
(the function you are trying to solve for) as ``func.func``. If, while trying
to solve the ODE, you find that you cannot, raise ``NotImplementedError``.
:py:meth:`~diofant.solvers.ode.dsolve` will catch this error with the ``all``
meta-hint, rather than causing the whole routine to fail.
Add a docstring to your function that describes the method employed. Like
with anything else in Diofant, you will need to add a doctest to the docstring,
in addition to real tests in ``test_ode.py``. Try to maintain consistency
with the other hint functions' docstrings. Add your method to the list at the
top of this docstring. Also, add your method to ``ode.rst`` in the
``docs/src`` directory, so that the Sphinx docs will pull its docstring into
the main Diofant documentation. Be sure to make the Sphinx documentation by
running ``make html`` from within the docs directory to verify that the
docstring formats correctly.
If your solution method involves integrating, use :py:meth:`Integral()
<diofant.integrals.integrals.Integral>` instead of
:py:meth:`~diofant.integrals.integrals.integrate`. This allows the user to bypass
hard/slow integration by using the ``_Integral`` variant of your hint. In
most cases, calling :py:meth:`diofant.core.basic.Basic.doit` will integrate your
solution. If this is not the case, you will need to write special code in
:py:meth:`~diofant.solvers.ode._handle_Integral`. Arbitrary constants should be
symbols named ``C1``, ``C2``, and so on. All solution methods should return
an equality instance. If you need an arbitrary number of arbitrary constants,
you can use ``constants = numbered_symbols(prefix='C', cls=Symbol, start=1)``.
If it is possible to solve for the dependent function in a general way, do so.
Otherwise, do as best as you can, but do not call solve in your
``ode_<hint>()`` function. :py:meth:`~diofant.solvers.ode.odesimp` will attempt
to solve the solution for you, so you do not need to do that. Lastly, if your
ODE has a common simplification that can be applied to your solutions, you can
add a special case in :py:meth:`~diofant.solvers.ode.odesimp` for it. For
example, solutions returned from the ``1st_homogeneous_coeff`` hints often
have many :py:meth:`~diofant.functions.elementary.exponential.log` terms, so
:py:meth:`~diofant.solvers.ode.odesimp` calls
:py:meth:`~diofant.simplify.simplify.logcombine` on them (it also helps to write
the arbitrary constant as ``log(C1)`` instead of ``C1`` in this case). Also
consider common ways that you can rearrange your solution to have
:py:meth:`~diofant.solvers.ode.constantsimp` take better advantage of it. It is
better to put simplification in :py:meth:`~diofant.solvers.ode.odesimp` than in
your method, because it can then be turned off with the simplify flag in
:py:meth:`~diofant.solvers.ode.dsolve`. If you have any extraneous
simplification in your function, be sure to only run it using ``if
match.get('simplify', True):``, especially if it can be slow or if it can
reduce the domain of the solution.
Finally, as with every contribution to Diofant, your method will need to be
tested. Add a test for each method in ``test_ode.py``. Follow the
conventions there, i.e., test the solver using ``dsolve(eq, f(x),
hint=your_hint)``, and also test the solution using
:py:meth:`~diofant.solvers.ode.checkodesol` (you can put these in a separate
tests and skip/XFAIL if it runs too slow/doesn't work). Be sure to call your
hint specifically in :py:meth:`~diofant.solvers.ode.dsolve`, that way the test
won't be broken simply by the introduction of another matching hint. If your
method works for higher order (>1) ODEs, you will need to run ``sol =
constant_renumber(sol, 'C', 1, order)`` for each solution, where ``order`` is
the order of the ODE. This is because ``constant_renumber`` renumbers the
arbitrary constants by printing order, which is platform dependent. Try to
test every corner case of your solver, including a range of orders if it is a
`n`\th order solver, but if your solver is slow, such as if it involves hard
integration, try to keep the test run time down.
Feel free to refactor existing hints to avoid duplicating code or creating
inconsistencies. If you can show that your method exactly duplicates an
existing method, including in the simplicity and speed of obtaining the
solutions, then you can remove the old, less general method. The existing
code is tested extensively in ``test_ode.py``, so if anything is broken, one
of those tests will surely fail.
"""
from collections import defaultdict
from itertools import islice
from ..core import (Add, AtomicExpr, Derivative, Dummy, E, Eq, Equality, Expr,
Function, I, Integer, Mul, Number, Pow, Subs, Symbol,
Tuple, Wild, diff, expand, expand_mul, factor_terms, nan,
oo, symbols, sympify, zoo)
from ..core.compatibility import is_sequence, iterable
from ..core.function import AppliedUndef, _mexpand
from ..core.multidimensional import vectorize
from ..functions import (atan2, conjugate, cos, exp, factorial, im, log, re,
sin, sqrt, tan)
from ..integrals import Integral, integrate
from ..logic.boolalg import BooleanAtom
from ..matrices import BlockDiagMatrix, Matrix, wronskian
from ..polys import Poly, PolynomialError, RootOf, lcm, terms_gcd
from ..polys.polyroots import roots_quartic
from ..polys.polytools import cancel, degree, div
from ..series import Order, series
from ..simplify import (collect, collect_const, cse, logcombine, posify,
powsimp, separatevars, simplify, trigsimp)
from ..utilities import default_sort_key, numbered_symbols, ordered, sift
from .deutils import _desolve, _preprocess, ode_order
from .pde import pdsolve
from .solvers import solve
#: This is a list of hints in the order that they should be preferred by
#: :py:meth:`~diofant.solvers.ode.classify_ode`. In general, hints earlier in the
#: list should produce simpler solutions than those later in the list (for
#: ODEs that fit both). For now, the order of this list is based on empirical
#: observations by the developers of Diofant.
#:
#: The hint used by :py:meth:`~diofant.solvers.ode.dsolve` for a specific ODE
#: can be overridden (see the docstring).
#:
#: In general, ``_Integral`` hints are grouped at the end of the list, unless
#: there is a method that returns an unevaluable integral most of the time
#: (which go near the end of the list anyway). ``default``, ``all``,
#: ``best``, and ``all_Integral`` meta-hints should not be included in this
#: list, but ``_best`` and ``_Integral`` hints should be included.
allhints = (
'separable',
'1st_exact',
'1st_linear',
'Bernoulli',
'Riccati_special_minus2',
'1st_homogeneous_coeff_best',
'1st_homogeneous_coeff_subs_indep_div_dep',
'1st_homogeneous_coeff_subs_dep_div_indep',
'almost_linear',
'linear_coefficients',
'separable_reduced',
'1st_power_series',
'lie_group',
'nth_linear_constant_coeff_homogeneous',
'nth_linear_euler_eq_homogeneous',
'nth_linear_constant_coeff_undetermined_coefficients',
'nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients',
'nth_linear_constant_coeff_variation_of_parameters',
'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters',
'Liouville',
'2nd_power_series_ordinary',
'2nd_power_series_regular',
'separable_Integral',
'1st_exact_Integral',
'1st_linear_Integral',
'Bernoulli_Integral',
'1st_homogeneous_coeff_subs_indep_div_dep_Integral',
'1st_homogeneous_coeff_subs_dep_div_indep_Integral',
'almost_linear_Integral',
'linear_coefficients_Integral',
'separable_reduced_Integral',
'nth_linear_constant_coeff_variation_of_parameters_Integral',
'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral',
'Liouville_Integral',
)
lie_heuristics = (
'abaco1_simple',
'abaco1_product',
'abaco2_similar',
'abaco2_unique_unknown',
'abaco2_unique_general',
'linear',
'function_sum',
'bivariate',
'chi'
)
def sub_func_doit(eq, func, new):
r"""
When replacing the func with something else, we usually want the
derivative evaluated, so this function helps in making that happen.
To keep subs from having to look through all derivatives, we mask them off
with dummy variables, do the func sub, and then replace masked-off
derivatives with their doit values.
Examples
========
>>> sub_func_doit(3*f(x).diff(x) - 1, f(x), x)
2
>>> sub_func_doit(x*f(x).diff(x) - f(x)**2 + f(x), f(x), 1/(x*(z + 1/x)))
x*(-1/(x**2*(z + 1/x)) + 1/(x**3*(z + 1/x)**2)) + 1/(x*(z + 1/x))
...- 1/(x**2*(z + 1/x)**2)
"""
reps = {}
repu = {}
for d in eq.atoms(Derivative):
u = Dummy('u')
repu[u] = d.subs({func: new}).doit()
reps[d] = u
return eq.subs(reps).subs({func: new.doit()}).subs(repu)
def get_numbered_constants(eq, num=1, start=1, prefix='C'):
"""
Returns a list of constants that do not occur
in eq already.
"""
if isinstance(eq, Expr):
eq = [eq]
elif not iterable(eq):
raise ValueError(f'Expected Expr or iterable but got {eq}')
atom_set = set().union(*[i.free_symbols for i in eq])
functions_set = set().union(*[i.atoms(Function) for i in eq])
if functions_set:
atom_set |= {Symbol(str(f.func)) for f in functions_set}
ncs = numbered_symbols(start=start, prefix=prefix, exclude=atom_set)
Cs = [next(ncs) for i in range(num)]
return Cs[0] if num == 1 else tuple(Cs)
def dsolve(eq, func=None, hint='default', simplify=True,
init=None, xi=None, eta=None, x0=0, n=6, **kwargs):
r"""
Solves any (supported) kind of ordinary differential equation and
system of ordinary differential equations.
**For single ordinary differential equation**
It is classified under this when number of equation in ``eq`` is one.
**Usage**
``dsolve(eq, f(x), hint)`` -> Solve ordinary differential equation
``eq`` for function ``f(x)``, using method ``hint``.
**Details**
``eq`` can be any supported ordinary differential equation (see the
:py:mod:`~diofant.solvers.ode` docstring for supported methods).
This can either be an :py:class:`~diofant.core.relational.Equality`,
or an expression, which is assumed to be equal to ``0``.
``f(x)`` is a function of one variable whose derivatives in that
variable make up the ordinary differential equation ``eq``. In
many cases it is not necessary to provide this; it will be
autodetected (and an error raised if it couldn't be detected).
``hint`` is the solving method that you want dsolve to use. Use
``classify_ode(eq, f(x))`` to get all of the possible hints for an
ODE. The default hint, ``default``, will use whatever hint is
returned first by :py:meth:`~diofant.solvers.ode.classify_ode`. See
Hints below for more options that you can use for hint.
``simplify`` enables simplification by
:py:meth:`~diofant.solvers.ode.odesimp`. See its docstring for more
information. Turn this off, for example, to disable solving of
solutions for ``func`` or simplification of arbitrary constants.
It will still integrate with this hint. Note that the solution may
contain more arbitrary constants than the order of the ODE with
this option enabled.
``xi`` and ``eta`` are the infinitesimal functions of an ordinary
differential equation. They are the infinitesimals of the Lie group
of point transformations for which the differential equation is
invariant. The user can specify values for the infinitesimals. If
nothing is specified, ``xi`` and ``eta`` are calculated using
:py:meth:`~diofant.solvers.ode.infinitesimals` with the help of various
heuristics.
``init`` is the set of initial/boundary conditions for the differential equation.
It should be given in the form of ``{f(x0): x1, f(x).diff(x).subs({x: x2}):
x3}`` and so on. For power series solutions, if no initial
conditions are specified ``f(0)`` is assumed to be ``C0`` and the power
series solution is calculated about 0.
``x0`` is the point about which the power series solution of a differential
equation is to be evaluated.
``n`` gives the exponent of the dependent variable up to which the power series
solution of a differential equation is to be evaluated.
**Hints**
Aside from the various solving methods, there are also some meta-hints
that you can pass to :py:meth:`~diofant.solvers.ode.dsolve`:
``default``:
This uses whatever hint is returned first by
:py:meth:`~diofant.solvers.ode.classify_ode`. This is the
default argument to :py:meth:`~diofant.solvers.ode.dsolve`.
``all``:
To make :py:meth:`~diofant.solvers.ode.dsolve` apply all
relevant classification hints, use ``dsolve(ODE, func,
hint="all")``. This will return a dictionary of
``hint:solution`` terms. If a hint causes dsolve to raise the
``NotImplementedError``, value of that hint's key will be the
exception object raised. The dictionary will also include
some special keys:
- ``order``: The order of the ODE. See also
:py:meth:`~diofant.solvers.deutils.ode_order` in
``deutils.py``.
- ``best``: The simplest hint; what would be returned by
``best`` below.
- ``best_hint``: The hint that would produce the solution
given by ``best``. If more than one hint produces the best
solution, the first one in the tuple returned by
:py:meth:`~diofant.solvers.ode.classify_ode` is chosen.
- ``default``: The solution that would be returned by default.
This is the one produced by the hint that appears first in
the tuple returned by
:py:meth:`~diofant.solvers.ode.classify_ode`.
``all_Integral``:
This is the same as ``all``, except if a hint also has a
corresponding ``_Integral`` hint, it only returns the
``_Integral`` hint. This is useful if ``all`` causes
:py:meth:`~diofant.solvers.ode.dsolve` to hang because of a
difficult or impossible integral. This meta-hint will also be
much faster than ``all``, because
:py:meth:`~diofant.core.expr.Expr.integrate` is an expensive
routine.
``best``:
To have :py:meth:`~diofant.solvers.ode.dsolve` try all methods
and return the simplest one. This takes into account whether
the solution is solvable in the function, whether it contains
any Integral classes (i.e. unevaluatable integrals), and
which one is the shortest in size.
See also the :py:meth:`~diofant.solvers.ode.classify_ode` docstring for
more info on hints, and the :py:mod:`~diofant.solvers.ode` docstring for
a list of all supported hints.
**Tips**
- See ``test_ode.py`` for many tests, which serves also as a set of
examples for how to use :py:meth:`~diofant.solvers.ode.dsolve`.
- :py:meth:`~diofant.solvers.ode.dsolve` always returns an
:py:class:`~diofant.core.relational.Equality` class (except for the
case when the hint is ``all`` or ``all_Integral``). If possible, it
solves the solution explicitly for the function being solved for.
Otherwise, it returns an implicit solution.
- Arbitrary constants are symbols named ``C1``, ``C2``, and so on.
- Because all solutions should be mathematically equivalent, some
hints may return the exact same result for an ODE. Often, though,
two different hints will return the same solution formatted
differently. The two should be equivalent. Also note that sometimes
the values of the arbitrary constants in two different solutions may
not be the same, because one constant may have "absorbed" other
constants into it.
- Do ``help(ode.ode_<hintname>)`` to get help more information on a
specific hint, where ``<hintname>`` is the name of a hint without
``_Integral``.
**For system of ordinary differential equations**
**Usage**
``dsolve(eq, func)`` -> Solve a system of ordinary differential
equations ``eq`` for ``func`` being list of functions including
`x(t)`, `y(t)`, `z(t)` where number of functions in the list depends
upon the number of equations provided in ``eq``.
**Details**
``eq`` can be any supported system of ordinary differential equations
This can either be an :py:class:`~diofant.core.relational.Equality`,
or an expression, which is assumed to be equal to ``0``.
``func`` holds ``x(t)`` and ``y(t)`` being functions of one variable which
together with some of their derivatives make up the system of ordinary
differential equation ``eq``. It is not necessary to provide this; it
will be autodetected (and an error raised if it couldn't be detected).
**Hints**
The hints are formed by parameters returned by classify_sysode, combining
them give hints name used later for forming method name.
Examples
========
>>> dsolve(Derivative(f(x), x, x) + 9*f(x), f(x))
Eq(f(x), C1*sin(3*x) + C2*cos(3*x))
>>> eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x)
>>> dsolve(eq, hint='1st_exact')
[Eq(f(x), -acos(C1/cos(x)) + 2*pi), Eq(f(x), acos(C1/cos(x)))]
>>> dsolve(eq, hint='almost_linear')
[Eq(f(x), -acos(C1/sqrt(-cos(x)**2)) + 2*pi), Eq(f(x), acos(C1/sqrt(-cos(x)**2)))]
>>> eq = (Eq(Derivative(f(t), t), 12*t*f(t) + 8*g(t)),
... Eq(Derivative(g(t), t), 21*f(t) + 7*t*g(t)))
>>> dsolve(eq)
[Eq(f(t), C1*x0(t) + C2*x0(t)*Integral(8*E**Integral(7*t, t)*E**Integral(12*t, t)/x0(t)**2, t)),
Eq(g(t), C1*y0(t) + C2*(E**Integral(7*t, t)*E**Integral(12*t, t)/x0(t) +
y0(t)*Integral(8*E**Integral(7*t, t)*E**Integral(12*t, t)/x0(t)**2, t)))]
>>> eq = (Eq(Derivative(f(t), t), f(t)*g(t)*sin(t)),
... Eq(Derivative(g(t), t), g(t)**2*sin(t)))
>>> dsolve(eq)
{Eq(f(t), -E**C1/(E**C1*C2 - cos(t))), Eq(g(t), -1/(C1 - cos(t)))}
"""
if iterable(eq):
match = classify_sysode(eq, func)
eq = match['eq']
order = match['order']
func = match['func']
t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
# keep highest order term coefficient positive
for i in range(len(eq)):
if eq[i].coeff(diff(func[i], t, ode_order(eq[i], func[i]))).is_negative:
eq[i] = -eq[i]
match['eq'] = eq
if len(set(order.values())) != 1:
raise ValueError('It solves only those systems of equations whose orders are equal')
match['order'] = list(order.values())[0]
if match['type_of_equation'] is None:
raise NotImplementedError
else:
if match['is_linear']:
if match['type_of_equation'] == 'type1' and match['order'] == 1:
solvefunc = globals()[f"sysode_linear_neq_order{match['order']}"]
elif match['no_of_equation'] <= 3:
solvefunc = globals()[f"sysode_linear_{match['no_of_equation']}eq_order{match['order']}"]
else:
raise NotImplementedError
else:
solvefunc = globals()[f"sysode_nonlinear_{match['no_of_equation']}eq_order{match['order']}"]
sols = solvefunc(match)
if init:
constants = Tuple(*sols).free_symbols - Tuple(*eq).free_symbols
solved_constants = solve_init(sols, func, constants, init)
return [sol.subs(solved_constants) for sol in sols]
return sols
else:
given_hint = hint # hint given by the user
# See the docstring of _desolve for more details.
hints = _desolve(eq, func=func,
hint=hint, simplify=True, xi=xi, eta=eta, type='ode', init=init,
x0=x0, n=n, **kwargs)
eq = hints.pop('eq', eq)
all_ = hints.pop('all', False)
if all_:
retdict = {}
gethints = classify_ode(eq, dict=True)
orderedhints = gethints['ordered_hints']
for hint in hints:
rv = _helper_simplify(eq, hint, hints[hint], simplify)
retdict[hint] = rv
func = hints[hint]['func']
retdict['best'] = min(list(retdict.values()), key=lambda x:
ode_sol_simplicity(x, func, trysolving=not simplify))
if given_hint == 'best':
return retdict['best']
for i in orderedhints:
if retdict['best'] == retdict.get(i, None):
retdict['best_hint'] = i
break
retdict['default'] = gethints['default']
retdict['order'] = gethints['order']
return retdict
else:
# The key 'hint' stores the hint needed to be solved for.
hint = hints['hint']
return _helper_simplify(eq, hint, hints, simplify, init=init)
def _helper_simplify(eq, hint, match, simplify=True, init=None, **kwargs):
r"""
Helper function of dsolve that calls the respective
:py:mod:`~diofant.solvers.ode` functions to solve for the ordinary
differential equations. This minimises the computation in calling
:py:meth:`~diofant.solvers.deutils._desolve` multiple times.
"""
r = match
if hint.endswith('_Integral'):
solvefunc = globals()['ode_' + hint[:-len('_Integral')]]
else:
solvefunc = globals()['ode_' + hint]
func = r['func']
order = r['order']
match = r[hint]
free = eq.free_symbols
def cons(s):
return s.free_symbols.difference(free)
if simplify:
# odesimp() will attempt to integrate, if necessary, apply constantsimp(),
# attempt to solve for func, and apply any other hint specific
# simplifications
sols = solvefunc(eq, func, order, match)
if isinstance(sols, Expr):
rv = odesimp(sols, func, order, cons(sols), hint)
else:
rv = [odesimp(s, func, order, cons(s), hint) for s in sols]
else:
# We still want to integrate (you can disable it separately with the hint)
match['simplify'] = False # Some hints can take advantage of this option
rv = _handle_Integral(solvefunc(eq, func, order, match),
func, order, hint)
if init and 'power_series' not in hint:
if isinstance(rv, Expr):
solved_constants = solve_init([rv], [r['func']], cons(rv), init)
rv = rv.subs(solved_constants)
elif iterable(rv):
rv1 = []
for s in rv:
solved_constants = solve_init([s], [r['func']], cons(s), init)
if solved_constants:
rv1.append(s.subs(solved_constants))
rv = rv1
else:
raise NotImplementedError
return rv
def solve_init(sols, funcs, constants, init):
"""
Solve for the constants given initial conditions
``sols`` is a list of solutions.
``funcs`` is a list of functions.
``constants`` is a list of constants.
``init`` is the set of initial/boundary conditions for the differential
equation. It should be given in the form of ``{f(x0): x1,
f(x).diff(x).subs({x: x2}): x3}`` and so on.
Returns a dictionary mapping constants to values.
``solution.subs(constants)`` will replace the constants in ``solution``.
Example
=======
>>> C1 = symbols('C1')
>>> sols = [Eq(f(x), C1*exp(x))]
>>> funcs = [f(x)]
>>> constants = [C1]
>>> init = {f(0): 2}
>>> solved_constants = solve_init(sols, funcs, constants, init)
>>> solved_constants
{C1: 2}
>>> sols[0].subs(solved_constants)
Eq(f(x), 2*E**x)
"""
# Assume init are of the form f(x0): value or Subs(diff(f(x), x, n), (x,
# x0)): value (currently checked by classify_ode). To solve, replace x
# with x0, f(x0) with value, then solve for constants. For f^(n)(x0),
# differentiate the solution n times, so that f^(n)(x) appears.
x = funcs[0].args[0]
diff_sols = []
subs_sols = []
diff_variables = set()
init = {k.doit() if isinstance(k, Subs) else k: v
for k, v in init.items()}
for funcarg, value in init.items():
if isinstance(funcarg, AppliedUndef):
x0 = funcarg.args[0]
matching_func = [f for f in funcs if f.func == funcarg.func][0]
S = sols
elif isinstance(funcarg, (Subs, Derivative)):
if isinstance(funcarg, Subs):
deriv = funcarg.expr
x0 = funcarg.point[0]
variables = funcarg.expr.variables
matching_func = deriv
else:
deriv = funcarg
x0 = funcarg.variables[0]
variables = (x,)*len(funcarg.variables)
matching_func = deriv.subs({x0: x})
if variables not in diff_variables:
for sol in sols:
if sol.has(deriv.expr.func):
diff_sols.append(Eq(sol.lhs.diff(*variables), sol.rhs.diff(*variables)))
diff_variables.add(variables)
S = diff_sols
else:
raise NotImplementedError('Unrecognized initial condition')
for sol in S:
if sol.has(matching_func):
sol2 = sol
sol2 = sol2.subs({x: x0})
sol2 = sol2.subs({funcarg: value})
subs_sols.append(sol2)
try:
solved_constants = solve(subs_sols, constants)
except NotImplementedError: # pragma: no cover
solved_constants = []
# XXX: We can't differentiate between the solution not existing because of
# invalid initial conditions, and not existing because solve is not smart
# enough. For now, we use NotImplementedError in this case.
if not solved_constants:
raise NotImplementedError("Couldn't solve for initial conditions")
if len(solved_constants) > 1: # pragma: no cover
raise NotImplementedError('Initial conditions produced too many solutions for constants')
return solved_constants[0]
def classify_ode(eq, func=None, dict=False, init=None, **kwargs):
r"""
Returns a tuple of possible :py:meth:`~diofant.solvers.ode.dsolve`
classifications for an ODE.
The tuple is ordered so that first item is the classification that
:py:meth:`~diofant.solvers.ode.dsolve` uses to solve the ODE by default. In
general, classifications at the near the beginning of the list will
produce better solutions faster than those near the end, thought there are
always exceptions. To make :py:meth:`~diofant.solvers.ode.dsolve` use a
different classification, use ``dsolve(ODE, func,
hint=<classification>)``. See also the
:py:meth:`~diofant.solvers.ode.dsolve` docstring for different meta-hints
you can use.
If ``dict`` is true, :py:meth:`~diofant.solvers.ode.classify_ode` will
return a dictionary of ``hint:match`` expression terms. This is intended
for internal use by :py:meth:`~diofant.solvers.ode.dsolve`. Note that
because dictionaries are ordered arbitrarily, this will most likely not be
in the same order as the tuple.
You can get help on different hints by executing
``help(ode.ode_hintname)``, where ``hintname`` is the name of the hint
without ``_Integral``.
See :py:data:`~diofant.solvers.ode.allhints` or the
:py:mod:`~diofant.solvers.ode` docstring for a list of all supported hints
that can be returned from :py:meth:`~diofant.solvers.ode.classify_ode`.
Notes
=====
These are remarks on hint names.
``_Integral``
If a classification has ``_Integral`` at the end, it will return the
expression with an unevaluated :py:class:`~diofant.integrals.integrals.Integral`
class in it. Note that a hint may do this anyway if
:py:meth:`~diofant.integrals.integrals.integrate` cannot do the integral,
though just using an ``_Integral`` will do so much faster. Indeed, an
``_Integral`` hint will always be faster than its corresponding hint
without ``_Integral`` because
:py:meth:`~diofant.integrals.integrals.integrate` is an expensive routine.
If :py:meth:`~diofant.solvers.ode.dsolve` hangs, it is probably because
:py:meth:`~diofant.core.expr.Expr.integrate` is hanging on a tough or
impossible integral. Try using an ``_Integral`` hint or
``all_Integral`` to get it return something.
Note that some hints do not have ``_Integral`` counterparts. This is
because :py:meth:`~diofant.integrals.integrals.integrate` is not used in solving
the ODE for those method. For example, `n`\th order linear homogeneous
ODEs with constant coefficients do not require integration to solve,
so there is no ``nth_linear_homogeneous_constant_coeff_Integrate``
hint. You can easily evaluate any unevaluated
:py:class:`~diofant.integrals.integrals.Integral`\s in an expression by doing
``expr.doit()``.
Ordinals
Some hints contain an ordinal such as ``1st_linear``. This is to help
differentiate them from other hints, as well as from other methods
that may not be implemented yet. If a hint has ``nth`` in it, such as
the ``nth_linear`` hints, this means that the method used to applies
to ODEs of any order.
``indep`` and ``dep``
Some hints contain the words ``indep`` or ``dep``. These reference
the independent variable and the dependent function, respectively. For
example, if an ODE is in terms of `f(x)`, then ``indep`` will refer to
`x` and ``dep`` will refer to `f`.
``subs``
If a hints has the word ``subs`` in it, it means the the ODE is solved
by substituting the expression given after the word ``subs`` for a
single dummy variable. This is usually in terms of ``indep`` and
``dep`` as above. The substituted expression will be written only in
characters allowed for names of Python objects, meaning operators will
be spelled out. For example, ``indep``/``dep`` will be written as
``indep_div_dep``.
``coeff``
The word ``coeff`` in a hint refers to the coefficients of something
in the ODE, usually of the derivative terms. See the docstring for
the individual methods for more info (``help(ode)``). This is
contrast to ``coefficients``, as in ``undetermined_coefficients``,
which refers to the common name of a method.
``_best``
Methods that have more than one fundamental way to solve will have a
hint for each sub-method and a ``_best`` meta-classification. This
will evaluate all hints and return the best, using the same
considerations as the normal ``best`` meta-hint.
Examples
========
>>> classify_ode(Eq(f(x).diff(x), 0), f(x))
('separable', '1st_linear', '1st_homogeneous_coeff_best',
'1st_homogeneous_coeff_subs_indep_div_dep',
'1st_homogeneous_coeff_subs_dep_div_indep',
'1st_power_series', 'lie_group',
'nth_linear_constant_coeff_homogeneous',
'separable_Integral', '1st_linear_Integral',
'1st_homogeneous_coeff_subs_indep_div_dep_Integral',
'1st_homogeneous_coeff_subs_dep_div_indep_Integral')
>>> classify_ode(f(x).diff(x, 2) + 3*f(x).diff(x) + 2*f(x) - 4)
('nth_linear_constant_coeff_undetermined_coefficients',
'nth_linear_constant_coeff_variation_of_parameters',
'nth_linear_constant_coeff_variation_of_parameters_Integral')
"""
init = sympify(init)
prep = kwargs.pop('prep', True)
if func and len(func.args) != 1:
raise ValueError('dsolve() and classify_ode() only '
f'work with functions of one variable, not {func}')
if prep or func is None:
eq, func_ = _preprocess(eq, func)
if func is None:
func = func_
x = func.args[0]
f = func.func
y = Dummy('y')
xi = kwargs.get('xi')
eta = kwargs.get('eta')
terms = kwargs.get('n')
if isinstance(eq, Equality):
if eq.rhs != 0:
return classify_ode(eq.lhs - eq.rhs, func, init=init, xi=xi,
n=terms, eta=eta, prep=False)
eq = eq.lhs
order = ode_order(eq, f(x))
# hint:matchdict or hint:(tuple of matchdicts)
# Also will contain "default":<default hint> and "order":order items.
matching_hints = {'order': order}
if not order:
if dict:
matching_hints['default'] = None
return matching_hints
else:
return ()
df = f(x).diff(x)
a = Wild('a', exclude=[f(x)])
b = Wild('b', exclude=[f(x)])
c = Wild('c', exclude=[f(x)])
d = Wild('d', exclude=[df, f(x).diff(x, 2)])
e = Wild('e', exclude=[df])
k = Wild('k', exclude=[df])
n = Wild('n', exclude=[f(x)])
c1 = Wild('c1', exclude=[x])
a2 = Wild('a2', exclude=[x, f(x), df])
b2 = Wild('b2', exclude=[x, f(x), df])
c2 = Wild('c2', exclude=[x, f(x), df])
d2 = Wild('d2', exclude=[x, f(x), df])
a3 = Wild('a3', exclude=[f(x), df, f(x).diff(x, 2)])
b3 = Wild('b3', exclude=[f(x), df, f(x).diff(x, 2)])
c3 = Wild('c3', exclude=[f(x), df, f(x).diff(x, 2)])
r3 = {'xi': xi, 'eta': eta} # Used for the lie_group hint
boundary = {} # Used to extract initial conditions
C1 = Symbol('C1')
eq = expand(eq)
# Preprocessing to get the initial conditions out
if init is not None:
init = {k.doit() if isinstance(k, Subs) else k: v
for k, v in init.items()}
for funcarg in init:
# Separating derivatives
if isinstance(funcarg, (Subs, Derivative)):
# f(x).diff(x).subs({x: 0}) is a Subs, but
# f(x).diff(x).subs({x: y}) is a Derivative
if isinstance(funcarg, Subs):
deriv = funcarg.expr
old = funcarg.variables[0]
new = funcarg.point[0]
else:
deriv = funcarg
# No information on this. Just assume it was x
old = x
new = funcarg.variables[0]
if (isinstance(deriv, Derivative) and
isinstance(deriv.args[0], AppliedUndef) and
deriv.args[0].func == f and len(deriv.args[0].args) == 1 and
old == x and not new.has(x) and
all(i == deriv.variables[0] for i in deriv.variables) and
not init[funcarg].has(f)):
dorder = ode_order(deriv, x)
temp = 'f' + str(dorder)
boundary.update({temp: new, temp + 'val': init[funcarg]})
else:
raise ValueError('Enter valid boundary conditions for Derivatives')
# Separating functions
elif isinstance(funcarg, AppliedUndef):
if (funcarg.func == f and len(funcarg.args) == 1 and
not funcarg.args[0].has(x) and not init[funcarg].has(f)):
boundary.update({'f0': funcarg.args[0], 'f0val': init[funcarg]})
else:
raise ValueError('Enter valid boundary conditions for Function')
else:
raise ValueError('Enter boundary conditions of the form init={f(point}: value, f(x).diff(x, order).subs({x: point}): value}')
# Precondition to try remove f(x) from highest order derivative
reduced_eq = None
if eq.is_Add:
deriv_coef = eq.coeff(f(x).diff(x, order))
if deriv_coef not in (1, 0):
r = deriv_coef.match(a*f(x)**c1)
if r and r[c1]:
den = f(x)**r[c1]
reduced_eq = Add(*[arg/den for arg in eq.args])
if not reduced_eq:
reduced_eq = eq
if order == 1: