Using sympy to create a matrix of symbolic coefficients
Zoho
I am brand new to python/sympy scripting and I have a project I have
written in Maxima but would like to port it to sympy (as it is far
more readable than the specific Maxima language) but I am not sure if
sympy yet has the funcialanity that I require.
What I would like to do is create an array with symbolic coefficients
and then use the array coeffiecents to create expressions. The
coeffiecents are solved for later. To give you and idea of what I mean
here is a small snippet of the Maxima code:
/* A Maxima program to solve
psi'' + lambda *sin(psi) = 0 , psi'(0)=psi'(1)=0,
int(sin(psi(s),s=0..1) = 0.
where psi = b*cos(n*pi*s) + b^3*cos^3(n*pi*s) + ...
lambda=1+k*b^2+k*b^4 + ... */
/* Setup the expressions for psi, lambda, x, y and solve for
the parameters */
setup(terms):=block([],
kill(k,a,b,d),
n:terms,
ratweight (b,1,c,1),
ratwtlvl: 2*(2*n+1),
/* Calculate expansion for psi and lambda */
C[0]:c,
Cpp[0]:-c,
k[0]:1,
for i thru n do (
C[i]:sum(a[i,j]*c^(2*j+1),j,1,i),
Cpp[i]:sum(a[i,j]*(-(2*j+1)^2*c^(2*j+1)+(2*j+1)*(2*j)*c^(2*j-1)),j,
1,i) ),
psi:sum(b^(2*i+1)*C[i], i,0,n),
lam:sum(b^(2*i)*k[i], i,0,n),
....
in the Maxima code you do not need to declare the array before. Here
'a' is an array and the expressions 'psi' and 'lam' are later solved
for the coefficients.
So is there a way to do this using sympy?
Zoho
post scriptum: Thank you for writing such great software!
Ondrej Certik
Hi Zoho,
thanks for your interest. I think you can do that in SymPy.
But I have problems understanding the Maxima code, could you please
write in python terms what it's doing? :)
For example:
> for i thru n do (
> C[i]:sum(a[i,j]*c^(2*j+1),j,1,i),
> Cpp[i]:sum(a[i,j]*(-(2*j+1)^2*c^(2*j+1)+(2*j+1)*(2*j)*c^(2*j-1)),j,
> 1,i) ),
This assignes some expressions into arrays C and Cpp, right?
> psi:sum(b^(2*i+1)*C[i], i,0,n),
> lam:sum(b^(2*i)*k[i], i,0,n),
This uses C to construct psi. Where is Cpp used? Isn't k[i]
uninitialized, except k[0]? What is in a[i,j],
if you didn't assign to it? But if you are just using arrays, symbolic
expressions, sums, substitution,
all of that work nice in sympy.
Ask here if you have any problems rewriting your code and please
report all bugs that you may find,
we'll try to fix them promptly.
Ondrej
Kirill Smelkov
Is this what you need?
========================================
#!/usr/bin/env python
from sympy import *
pprint_try_use_unicode()
c = Symbol('c')
b = Symbol('b')
def irange(a,b):
"""inclusive range -- range(a,b+1)"""
return range(a,b+1)
def setup_terms(n):
a = {}
C = [None]*(n+1)
C[0] = c
for i in irange(1,n):
for j in irange(1,i):
a[i,j] = Symbol('a_%i%i' % (i,j))
C[i] = sum(a[i,j] * c**(2*j+1) for j in irange(1,i))
psi = sum(b**(2*i+1) * C[i] for i in irange(0,n))
return psi
if __name__ == '__main__':
psi = setup_terms(3)
print 'PSI:'
pprint(psi)
========================================
On my host the output is:
kirr@evo:~/src/sympy/sympy$ ./psi.py
PSI:
5 ⎛ 3 5⎞ 7 ⎛ 3 5 7⎞ 3 3
b*c + b *⎝a₂₁*c + a₂₂*c ⎠ + b *⎝a₃₁*c + a₃₂*c + a₃₃*c ⎠ + a₁₁*b *c
> post scriptum: Thank you for writing such great software!
Thanks for your interest to SymPy too!
--
Всего хорошего, Кирилл.
http://landau.phys.spbu.ru/~kirr/aiv/
Zoho
what the code actually does, I should have posted it all (and I will
following this brief explanation). I am using a modified perturbation
expansion technique in order to solve a particular 2nd order nonlinear
ODE. The equation describes how a beam responds to a force applied to
its endpoints where the line of action of the force, initially
directed along the axis of the beam, does not change.
Here 's' describes a position along the beam, 'psi(s)' is the angle
that the beam makes with the line of action of the force and 'lambda'
is proportional to the applied force. I assume that psi can be
described as: psi = b*cos(n*pi*s) + b**3*a_11*cos**3(n*pi*s) +
b**5*(a_21*cos**3(n*pi*s) + a_22*cos**5(n*pi*s)) + ... and the
parameter 'lambda' can be described as lam = n*pi**2*(1 + k_1*b**2 +
k_2*b**4 + ...)
So if we can accept that, then the perturbation code takes as its
argument the number of terms wanted in the expansion and then
calulates the values of k and a. After k and a have been calculated
we can pass a value for 'b' (which is arguably analgous to the initate
angle of the beam) and then convert the solution for 'psi' into
functions for x and y. Finally plot the result.
Now, in the orginal code I wrote I used a few built in functions
available in Maxima which I cannot find the proper syntax for in
sympy. I need to be able to simply an expression around a specfic
variable and then extract the coefficeints for said variable. Here is
an example of what I mean:
/* Solve for parameters */
for j:0 thru i-1 do (
map(rhs, solve(ratcoef(A,c,2*(i-j)+1),a[i,i-j])),
a[i,i-j]:first(ev(%%)) ),
map(rhs, solve(ratcoef(A,c,1),k[i])),
k[i]:first(ev(%%)) ),
here we can see in the declaration for A we use Cpp (where the p
stands for prime so Cpp is the second derivative of C wrt 's') and add
it to the specific term in lam*sin*psi. That allows us to solve for k
and a. So ratcoef(lamsinpsi,b,2*i+1) takes the expression lamsinpsi
and organizes it around powers of b, and then passes out the
coefficient of the power of b requested.
Can anyone point me in a quasi-proper direction for preforming the
same calculations in sympy? Thanks agian for all your help!
Zoho
/* Setup the expressions for psi, lambda, x, y and solve for
the parameters */
setup(terms):=block([],
kill(k,a,b,d),
n:terms,
ratweight (b,1,c,1),
ratwtlvl: 2*(2*n+1),
/* Calculate expansion for psi and lambda */
C[0]:c,
Cpp[0]:-c,
k[0]:1,
for i thru n do (
C[i]:sum(a[i,j]*c^(2*j+1),j,1,i),
Cpp[i]:sum(a[i,j]*(-(2*j+1)^2*c^(2*j+1)+(2*j+1)*(2*j)*c^(2*j-1)),j,
1,i) ),
psi:sum(b^(2*i+1)*C[i], i,0,n),
lam:sum(b^(2*i)*k[i], i,0,n),
sinpsi:sum((-1)^i/(2*i+1)!*psi^(2*i+1),i,0,n),
/* Calculate lam*sin(psi) */
lamsinpsi:ratexpand(lam*sinpsi),
/* Solve for parameters */
for j:0 thru i-1 do (
map(rhs, solve(ratcoef(A,c,2*(i-j)+1),a[i,i-j])),
a[i,i-j]:first(ev(%%)) ),
map(rhs, solve(ratcoef(A,c,1),k[i])),
k[i]:first(ev(%%)) ),
/* Calculate cos(psi) = 1 - psi^2/2! + ... */
ratwtlvl:4*n,
cospsi:sum((-1)^i/(2*i)!*psi^(2*i),i,1,n)+1,
/* Calculate the integrals int_0^s cos(psi) and int_0^s sin(psi) */
cospsi:ev(ratexpand(cospsi)),
ratwtlvl:2*(2*n+1),
sinpsi:ev(ratexpand(sinpsi)),
ratwtlvl:false,
x(s):=sum(ratcoef(cospsi,c,2*i)*integrate(cos(m*%pi*s)^(2*i),s),i,
0,n),
y(s):=sum(ratcoef(sinpsi,c,2*i+1)*integrate(cos(m*%pi*s)^(2*i+1),s),i,
0,n) ) $
/* Plot a graph of the beam in the xy-plane */
beam(b,m):=block([s],
plot2d([parametric,x(s),y(s),[s,0,1]],[nticks,200]) ) $
Ondrej Certik
Hi Zoho,
first I am sorry for my later reply, I forgot about this.
Unfortunately, I don't have time right now to learn the syntax of
maxima. Could you please post here some examples of usage of the above
code? I.e. which expression you plug in and what you get? It will be
then easy for me to help you do the same in SymPy.
Thanks a lot,
Ondrej
Zoho Vignochi
Actually I have altered the code slightly.I believe the only replacement
routine I need now is for ratcoef (in Maxima) or coeff and collect (in
Maple). I guess it is easiest to explain what the maple commands do and
then see if that is possible in sympy.
syntax coeff(p,x,n)
The coeff command takes a polynomial p and picks out the coefficient of
the term containing the variable x to the given power n.
An example for coeff:
f := (a+1)*x + (a+2)*x^2 + a;
-> f := (a+1)x + (a+2)x^2 + a
coeff(f,x,2);
-> a+2
coeff(f,x,0)
-> a
And an example for collect:
syntax collect(p,x)
The collect command is used to collect like terms of the variable x in
the polynomial p.
f := (a+1)*x + a^2*x + a;
-> f := (a+1)x + a^2*x + a
collect(f,x);
-> (a+1 + a^2)x + a
collect(f,a);
-> a^2*x + (x+1)a + x
I hope that helps. I really appreciate your help. I have read through
the documentation on the website but didn't find answers there. Is there
another source of material available to read through? Perhaps these are
really simple questions and I just didn't find the right things to read!
Thanks again,
Zoho
Ondrej Certik
You can use this:
In [1]: a = Symbol("a")
In [2]: f = (a+1)*x + (a+2)*x**2 + a
In [3]: f = Polynomial(f, var=[x])
In [4]: f.nth_coeff(2)
Out[4]: 2 + a
In [5]: f.nth_coeff(0)
Out[5]: a
But using the method you suggested is a more common interface, so I
implemented that:
http://code.google.com/p/sympy/issues/detail?id=691
But it still needs some polishing.
>
> And an example for collect:
> syntax collect(p,x)
>
> The collect command is used to collect like terms of the variable x in
> the polynomial p.
>
> f := (a+1)*x + a^2*x + a;
> -> f := (a+1)x + a^2*x + a
> collect(f,x);
> -> (a+1 + a^2)x + a
> collect(f,a);
> -> a^2*x + (x+1)a + x
This is implemented using the same syntax:
In [1]: a = Symbol("a")
In [2]: f = (a+1)*x + a**2*x+a
In [3]: collect(f, x)
Out[3]:
⎛ 2⎞
a + x*⎝1 + a + a ⎠
In [4]: collect(f, a)
Out[4]:
2
a + x*a + x*(1 + a)
>
> I hope that helps. I really appreciate your help. I have read through
> the documentation on the website but didn't find answers there. Is there
> another source of material available to read through? Perhaps these are
> really simple questions and I just didn't find the right things to read!
The best thing is to read through sources (all methods should have
understandable docstrings) and especially the tests, because all
features have tests in sympy, so just by reading those you get the
idea what you can do. Another good approach is to start isympy and hit
"coll" + TAB, this completes to "collect", then you type "collect?"
and it will give you a docstring.
Ondrej
Ondrej Certik
> > f := (a+1)*x + (a+2)*x^2 + a;
> > -> f := (a+1)x + (a+2)x^2 + a
> > coeff(f,x,2);
> > -> a+2
> > coeff(f,x,0)
> > -> a
.coeff() was just pushed into our main hg repo, so update your sympy and do:
In [1]: a = Symbol("a")
In [2]: f = (a+1)*x + (a+2)*x**2 + a
In [3]: f = Polynomial(f)
In [4]: f.coeff(x, 2)
Out[4]: 2 + a
In [5]: f.coeff(x, 0)
Out[5]: a
As always, documentation can be found using for example:
In [6]: f.coeff?
Type: instancemethod
Base Class: <type 'instancemethod'>
String Form:
<bound method Polynomial.coeff of Polynomial( 2 2
a + x + 2*x + a*x + a*x , ((1, 1, 2), (1, 1, 1), (2, 0, 2), (1, 1,
0), (1, 0, 1)), [a, x], 'grevlex')>
Namespace: Interactive
File: /home/ondra/sympy/sympy/polynomials/base.py
Definition: f.coeff(self, x, n)
Docstring:
Returns the coefficient at x**n
Example:
>>> a, x = symbols("ax")
>>> f = (a+1)*x + (a+2)*x**2 + a
>>> Polynomial(f).coeff(x, 2)
2 + a
>>> Polynomial(f).coeff(a, 1)
1 + x + x**2
The messy output of repr(f) is due to:
http://code.google.com/p/sympy/issues/detail?id=659
Ondrej
Zoho Vignochi
I do feel embarrassed that I asked about collect but it was already
there! Sorry for the unnecessary noise.
I will update sympy now to get .coeff and thank you so much for all your
help!
Zoho
Ondrej Certik
> there! Sorry for the unnecessary noise.
No problem, that's what the mailinglist is for. Also I didn't know it
either off top of my head. What I do is that I write the expression
into a variable in isympy and then
t.<TAB>
this prints all the available method names. If I didn't find the one I
want, I do:
dir()
this prints all functions, that are imported, so I do
In [2]: c<TAB>
cache chebyshevu_root combsimp cos
callable chr compile cosh
cancel class complex cot
casoratian classmethod concrete coth
cat clear conjugate count_real_roots
ceiling cmp continue cp
chebyshevt coerce convex_hull credits
chebyshevt_root collect copyright
chebyshevu combinatorial core
and that's it. (Then collect? and it will print a nice help with
examples of usage).
> I will update sympy now to get .coeff and thank you so much for all your
> help!
Thanks for your feedback! Please write when you have more problems/questions.
Ondrej