2008/5/27 Friedrich Hagedorn <[EMAIL PROTECTED]>:
> On Tue, May 27, 2008 at 12:07:33PM +0200, Friedrich Hagedorn wrote:
>>
>> On Tue, May 27, 2008 at 11:42:21AM +0200, Ondrej Certik wrote:
>> > In [4]: g_dd
>> > Out[4]:
>> > ⎡ ⎛ 2 2⎞ ⎤
>> > ⎢-1 + ω*⎝x + y ⎠ 2*ω*y -2*ω*x 0⎥
>> > ⎢ ⎥
>> > ⎢ 2*ω*y 1 0 0⎥
>> > ⎢ ⎥
>> > ⎢ -2*ω*x 0 1 0⎥
>> > ⎢ ⎥
>> > ⎣ 0 0 0 1⎦
> [...]
>> > In [3]: g_uu = g_dd.inv()
> [...]
>> I guess the matrix inversion is done by the gauss algorithm. But my
>> experience is that the gauss isnt so good for symbolic calculation.
>>
>> Try to calc the matrix inverse with the determiants
>>
>> M^-1 = 1/det(M) * (M_adj).T
>>
>> with
>>
>> M_adj := (-1)^(i+j) * det(M_ij)
>>
>> M_ij := Matrix M without row i and without column j
>>
>> and the evaluations of the determinat of your matrix is easy. This is
>> also (shortly) discribed in
>
> So here is my quick solution with the attached file inversion.py:
>
> In [1]: M = Matrix([(z,y,-x,0),(y, 1, 0, 0), (-x, 0, 1, 0), (0, 0, 0, 1)])
>
> In [2]: M
> Out[2]:
> ⎡ z y -x 0⎤
> ⎢ y 1 0 0⎥
> ⎢-x 0 1 0⎥
> ⎣ 0 0 0 1⎦
>
> In [3]: run inverse.py
>
> In [4]: adj(M).T
> Out[4]:
> ⎡ ⎤
> ⎢ 1 -y x 0⎥
> ⎢ 2 ⎥
> ⎢ -y z - x -x*y 0⎥
> ⎢ 2 ⎥
> ⎢ x -x*y z - y 0⎥
> ⎢ 2 2⎥
> ⎣ 0 0 0 z - x - y ⎦
>
> In [5]: 1/M.det()
> Out[5]:
> 1
> ───────────
> 2 2
> z - x - y
>
> If you multiply Out[4] * Out[5] you get the inversion matrix and can test
>
> In [15]: (M*adj(M).T)*(1/M.det())
> Out[15]:
> ⎡1 0 0 0⎤
> ⎢0 1 0 0⎥
> ⎢0 0 1 0⎥
> ⎣0 0 0 1⎦
>
> I hope that helps.
Wow, indeed, that helps a lot!
With my original problem, I am getting this:
In [18]: e = adj(M).T * (1/M.det())
In [19]: one = M*e
In [20]: one[0,0]
Out[20]:
2
1 ω*x
- ──────────────────────────────────── + ────────────────────────────────────
2 2 2 2 2 2 2 2 2 2 2 2
-1 + ω*x + ω*y - 4*ω *x - 4*ω *y -1 + ω*x + ω*y - 4*ω *x - 4*ω *y
2 2 2
ω*y 4*ω *x
+ ──────────────────────────────────── - ────────────────────────────────────
2 2 2 2 2 2 2 2 2 2 2 2
-1 + ω*x + ω*y - 4*ω *x - 4*ω *y -1 + ω*x + ω*y - 4*ω *x - 4*ω *y
2 2
4*ω *y
- ────────────────────────────────────
2 2 2 2 2 2
-1 + ω*x + ω*y - 4*ω *x - 4*ω *y
In [21]: simplify(one[0,0])
Out[21]: 1
So first the result looks much nicer and second simplify can handle
that. Very cool, I created a new issue for your code to get in:
http://code.google.com/p/sympy/issues/detail?id=866
Ondrej
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sympy?hl=en
-~----------~----~----~----~------~----~------~--~---
