how reliable is sympy in solving systems of differential equations

[krastano...@gmail.com]

Is there a straightforward way to solve this system in sympy:

simple linear system of linear equations with initial conditions

[f_1(t) + Derivative(f_0(t), t), -f_0(t) + Derivative(f_1(t), t),
f_0(0) - 1, f_1(0) - 1]

dsolve does not permit this directly.

[Aaron Meurer]

Not yet, I'm afraid, though most of the work that would be needed to
do it is already there (the exception is for non-diagonalizable
systems, for which the methods are not fully implemented in the
matrices yet).

By the way, it never occurred to me to allow initial conditions to
just be entered as equations. My idea was to enter them as a
dictionary, like dsolve([f(x).diff(x) - f(x)], ics={f(0)=1}). Do you
think that doing it your way would be better? I suppose it would be
possible to be fancy and define coupled initial conditions. See issue
1621.

Aaron Meurer

On May 26, 2012, at 5:01 PM, "krastano...@gmail.com"

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
>

[krastano...@gmail.com]

On 28 May 2012 00:34, Aaron Meurer <asme...@gmail.com> wrote:
> Not yet, I'm afraid, though most of the work that would be needed to
> do it is already there (the exception is for non-diagonalizable
> systems, for which the methods are not fully implemented in the
> matrices yet).
>

If I want to do it "by hand" what are the steps (I am asking in order
not to miss something that is already implemented)?

> By the way, it never occurred to me to allow initial conditions to
> just be entered as equations.  My idea was to enter them as a
> dictionary, like dsolve([f(x).diff(x) - f(x)], ics={f(0)=1}).  Do you
> think that doing it your way would be better?  I suppose it would be
> possible to be fancy and define coupled initial conditions. See issue
> 1621.

It is just a habit I have from Mathematica. I suppose that there are
some really weird cases where it may be advantageous, but for now I
prefer your proposal about kwargs.

[Aaron Meurer]

On May 28, 2012, at 4:26 AM, "krastano...@gmail.com"

<krastano...@gmail.com> wrote:

> On 28 May 2012 00:34, Aaron Meurer <asme...@gmail.com> wrote:
>> Not yet, I'm afraid, though most of the work that would be needed to
>> do it is already there (the exception is for non-diagonalizable
>> systems, for which the methods are not fully implemented in the
>> matrices yet).
>>
>
> If I want to do it "by hand" what are the steps (I am asking in order
> not to miss something that is already implemented)?

1. Convert the system to a matrix.
2. Diagonalize the matrix.
3. Solve the diagonal system.
4. Find the initial conditions in the diagonal basis.
5. Solve for the arbitrary constants.
6. Put the solution back in the original basis.

Currently only 2 and 3 are implemented (3 is of course trivial). Note
that if you want you can also do 1, 2, 3, 6, 5.

If you're interested in extending dsolve(), there's an extensive
docstring at the top of ode.py explaining everything you need to do.
The most difficult part is actually not what you expect: matching the
ODE and pulling out the coefficients in a general way.

>
>> By the way, it never occurred to me to allow initial conditions to
>> just be entered as equations. My idea was to enter them as a
>> dictionary, like dsolve([f(x).diff(x) - f(x)], ics={f(0)=1}). Do you
>> think that doing it your way would be better? I suppose it would be
>> possible to be fancy and define coupled initial conditions. See issue
>> 1621.
>
> It is just a habit I have from Mathematica. I suppose that there are
> some really weird cases where it may be advantageous, but for now I
> prefer your proposal about kwargs.

Well now I'm starting to like your way. Your way allows much more
general things, though it also puts a lot more parsing burden on the
implementation.

Aaron Meurer

[krastano...@gmail.com]

I have one more question about this part:

> 1. Convert the system to a matrix.

Do some other parts of sympy do this in a way that can catch nonlinear
equations that can trivially be transformed into linear ones. For
instance log(f'(t))=log(t) -> f'(t)=t.

Obviously sometimes the transformations are valid only over a certain
domain and probably there are other issues but maybe `solve()` already
does some of this magic.

[Aaron Meurer]

Not that I know of, but that doesn't mean that such a function doesn't
exist. solve() sounds like a good place to look.

You are starting to see why I said that the matching step is the
hardest. When working with the ODE module, you soon start to see the
limitations of pattern matching in SymPy.

Aaron Meurer

On May 28, 2012, at 1:47 PM, "krastano...@gmail.com"

[Chris Smith]

>> Do some other parts of sympy do this in a way that can catch nonlinear
>> equations that can trivially be transformed into linear ones. For
>> instance log(f'(t))=log(t) -> f'(t)=t.
>>

Solve solves this:

>>> solve(log(f(t).diff(t))-log(t), f(t).diff(t))
[t]

/c

[krastano...@gmail.com]

I encountered a problem when implementing this.

At the moment the constant of integration are Symbols and not Dummies.
Thus there is a clash between constants from different equations. How
appropriate is to make these constant Dummies?

[krastano...@gmail.com]

If I simply substitute Dummy for Symbol there are test failures
because the tests depend on Symbol('C1')==Symbol('C2').

On 30 May 2012 12:30, krastano...@gmail.com

[krastano...@gmail.com]

How appropriate it is to create a very simple IntegrationConstant
subclass of Symbol?

On 30 May 2012 12:39, krastano...@gmail.com

[krastano...@gmail.com]

On 30 May 2012 12:41, krastano...@gmail.com

<krastano...@gmail.com> wrote:
> How appropriate it is to create a very simple IntegrationConstant
> subclass of Symbol?

To be clearer: IntegrationConstant will have two internal counters:
equation_counter, to say whether the integration constants are
concerned with the same equation and visual counter that will be
printed (the second one in not necessary. it can be just a part of the
name as it is now)

[krastano...@gmail.com]

So here is a first try at solving this:

https://github.com/sympy/sympy/pull/1322

Comments will be very much appreciated.

[Aaron Meurer]

There is an issue about this:
http://code.google.com/p/sympy/issues/detail?id=1739

Aaron Meurer