collecting on similar series expansion
Hi , I'm just trying out ginac for the first time. I have an expression with sumands of this type (in latex format): 1-\frac{1}{3} \frac{ p D[0](F)(q) D[0,0](F)(q) \nu \epsilon^{2} dt^{5} \delta}{m^{3}}-\frac{1}{9} \frac{ p D[0](F)(q) D[0,0](F)(q)^{2} \alpha^{2} \nu^{2} \epsilon^{2} dt^{9} \F(q) \delta \gamma}{m^{5}}-\frac{1}{18} \frac{ p^{2} D[0,0](F)(q) D[0,0,0](F)(q) \beta \alpha^{2} \epsilon^{2} dt^{8} \F(q) \delta \gamma}{m^{5}}+\frac{1}{2} \frac{ D[0](F)(q) D[0,0](F)(q) \alpha \nu \epsilon dt^{6} \F(q)}{m^{3}} .... so i want to collect terms that have the same coefficients of the form D[0..i times] (F)(q) * D[ 0.. j times] (F)(q) .... but the polynomial collect won't work. Initially i thought about substituting each D[ ] (F)(q) with a different symbol, then collect on the multivariate polynomials of these symbols and then finally substitute back the D[] expressions, however, I don't know how the collect call is supposed to receive the symbols that it must collect, the code looks like this: ex jacobian = diff( q + dq_new*dt , q )*diff( p + dp_new*dt , p ) - diff( q + dq_new*dt , p )*diff( p + dp_new*dt , q ); cout << " the jacobian of the integrator is " << expand(jacobian) << endl; exmap diff_map , inv_diff_map; symbol foo[10]; diff_map[ F(q) ] = foo[0]; inv_diff_map[ foo[0] ] = F(q); diff_map[ F(q).diff(q) ] = foo[1]; inv_diff_map[ foo[1] ] = F(q).diff(q); diff_map[ F(q).diff(q,2) ] = foo[2]; inv_diff_map[ foo[2] ] = F(q).diff(q,2); diff_map[ F(q).diff(q,3) ] = foo[3]; inv_diff_map[ foo[3] ] = F(q).diff(q,3); diff_map[ F(q).diff(q,4) ] = foo[4]; inv_diff_map[ foo[4] ] = F(q).diff(q,4); ex rord = collect( jacobian.subs( diff_map ) , { foo[0] , foo[1] , foo[2] , foo[3] , foo[4] } ); // i know this is an error in C++, but this is the idea rord = rord.subs( inv_diff_map ); cout << "reordered we obtain " << rord << endl; any thoughts on this? __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas
Hello, On Sun, Jan 28, 2007 at 09:31:18AM -0300, Charlls Quarra wrote:
I'm just trying out ginac for the first time. I have an expression with sumands of this type (in latex format): [snipped the expression]
so i want to collect terms that have the same coefficients of the form D[0..i times] (F)(q) * D[ 0.. j times] (F)(q) .... but the polynomial collect won't work.
Could you please post the complete example which illustrates this?
Initially i thought about substituting each D[ ] (F)(q) with a different symbol, then collect on the multivariate polynomials of these symbols and then finally substitute back the D[] expressions, however, I don't know how the collect call is supposed to receive the symbols that it must collect,
Reading the manual _does_ help: "Another useful representation of multivariate polynomials is as a univariate polynomial in one of the variables with the coefficients being polynomials in the remaining variables. The method `collect()' accomplishes this task: ex ex::collect(const ex & s, bool distributed = false); The first argument to `collect()' can also be a list of objects in which case the result is either a recursively collected polynomial, or a polynomial in a distributed form with terms like c*x1^e1*...*xn^en, as specified by the `distributed' flag. Note that the original polynomial needs to be in expanded form (for the variables concerned) in order for `collect()' to be able to find the coefficients properly."
the code looks like this:
ex jacobian = diff( q + dq_new*dt , q )*diff( p + dp_new*dt , p ) - diff( q + dq_new*dt , p )*diff( p + dp_new*dt , q ); cout << " the jacobian of the integrator is " << expand(jacobian) << endl; exmap diff_map , inv_diff_map; symbol foo[10]; diff_map[ F(q) ] = foo[0]; inv_diff_map[ foo[0] ] = F(q); diff_map[ F(q).diff(q) ] = foo[1]; inv_diff_map[ foo[1] ] = F(q).diff(q); diff_map[ F(q).diff(q,2) ] = foo[2]; inv_diff_map[ foo[2] ] = F(q).diff(q,2); diff_map[ F(q).diff(q,3) ] = foo[3]; inv_diff_map[ foo[3] ] = F(q).diff(q,3); diff_map[ F(q).diff(q,4) ] = foo[4]; inv_diff_map[ foo[4] ] = F(q).diff(q,4); ex rord = collect( jacobian.subs( diff_map ) , { foo[0] , foo[1] , foo[2] , foo[3] , foo[4] } ); // i know this is an error in C++, but this is the idea rord = rord.subs( inv_diff_map ); cout << "reordered we obtain " << rord << endl;
any thoughts on this?
I tnink this should do the job: ex rord = jacobian.expand().collect(lst(F(q), F(q).diff(q), F(q).diff(q, 2), F(q).diff(q, 3), F(q).diff(q, 4)), true /* treat the expression as Z[x,y,...] (as opposed to Z[x][y]...) */ ); Best regards, Alexei -- All science is either physics or stamp collecting.
--- Sheplyakov Alexei <varg@theor.jinr.ru> escribió:
Could you please post the complete example which illustrates this?
I tnink this should do the job:
ex rord = jacobian.expand().collect(lst(F(q), F(q).diff(q), F(q).diff(q, 2), F(q).diff(q, 3), F(q).diff(q, 4)), true /* treat the expression as Z[x,y,...] (as opposed to Z[x][y]...) */ );
now seeing that example (and trying it out) i realize what was my error; i was calling jacobian.expand() and assumed that was enough to change the jacobian expression, but when i called collect later on the (still unexpanded) jacobian ex it wasn't doing the job at all. thanks this is what i was doing #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; DECLARE_FUNCTION_1P(F) static ex F_eval(const ex & x) { return F(x).hold(); } static ex F_evalf(const ex & x) { return F(x).hold(); } REGISTER_FUNCTION(F, eval_func(F_eval). evalf_func(F_evalf). latex_name("\\F")); int main() { cout << latex; symbol x("x"), y("y"); // ex F; ex f = pow( x , 5 )*sin( x )+ F(x); cout << " f is " << f << endl; cout << " D(f) is " << f.diff(x) << endl; cout << " D(f,2) is " << f.diff(x,2) << endl; cout << " integrator " << endl; symbol q("q") , p("p"); symbol m("m"), dt("dt"); symbol alpha("alpha") , beta("beta") , gamma("gamma") , delta("delta"), epsilon("epsilon") , nu("nu"); symbol dq_iota("dq_iota"); symbol dq_lambda("dq_lambda"); ex dq = p/m - F(q)*dt/m - F(q).diff(q)*dq_iota*pow(dt,2)/(2*m) - gamma*(F(q).diff(q,2)*pow(dq_iota,2) + F(q).diff(q)* dq_iota.diff(q) )*pow(dt,3)/(factorial(3) * m); ex dp = -F(q) -F(q).diff(q)*dq_lambda*dt - delta*(F(q).diff(q,2)*pow(dq_lambda,2) + F(q).diff(q)*dq_lambda.diff(q) )*pow( dt , 2 )/factorial(3); cout << "dq is " << dq << endl; cout << "dp is " << dp << endl; exmap s_map; s_map[ dq_iota ] = alpha*(p/m - beta*F(q)*dt/m); s_map[ dq_lambda ] = epsilon*(p/m - nu*F(q)*dt/m); ex dq_new = dq; ex dp_new = dp; dq_new = dq_new.subs( s_map ); dp_new = dp_new.subs( s_map ); cout << "now dq is " << dq_new << endl; cout << "now dp is " << dp_new << endl; ex jacobian = diff( q + dq_new*dt , q )*diff( p + dp_new*dt , p ) - diff( q + dq_new*dt , p )*diff( p + dp_new*dt , q ); cout << " the jacobian of the integrator is " << expand(jacobian) << endl; exmap diff_map , inv_diff_map; symbol foo0 , foo1 , foo2 , foo3 , foo4; lst params; params = foo0 , foo1 , foo2 , foo3 , foo4; diff_map[ F(q) ] = foo0; inv_diff_map[ foo0 ] = F(q); diff_map[ F(q).diff(q) ] = foo1; inv_diff_map[ foo1 ] = F(q).diff(q); diff_map[ F(q).diff(q,2) ] = foo2; inv_diff_map[ foo2 ] = F(q).diff(q,2); diff_map[ F(q).diff(q,3) ] = foo3; inv_diff_map[ foo3 ] = F(q).diff(q,3); diff_map[ F(q).diff(q,4) ] = foo4; inv_diff_map[ foo4 ] = F(q).diff(q,4); ex rord = collect( jacobian.subs( diff_map ) , params ); rord = rord.subs( inv_diff_map ); cout << "reordered we obtain " << rord << endl; ex eqns = jacobian == 1; ex vars = alpha; cout << lsolve( eqns , vars ) << endl; }; __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas
--- Charlls Quarra <charlls_quarra@yahoo.com.ar> escribió:
now seeing that example (and trying it out) i realize what was my error; i was calling jacobian.expand() and assumed that was enough to change the jacobian expression, but when i called collect later on the (still unexpanded) jacobian ex it wasn't doing the job at all.
The other problem i had with this program is that there wasnt a builtin way to get coefficients from multivariate polynomials, and i thought that i could do that in a simple manner with collect and subs, but i got this map_function which does what i want (run the program): the degree_list function is really meant to be a general all-in-once degree/coeff. i think that a function like this could be more optimal than to check full ranges between ldegree to degree, just a list with the available degrees.. just a thought #include <iostream> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; lst degree_list(const ex& poly_ex ,const ex& var) { lst result = lst(); ex cf; for (int i=poly_ex.ldegree( var ) ; i <= poly_ex.degree( var ) ; i++) { cf = poly_ex.coeff( var , i ); if (cf != 0) { result.append( lst( i , cf ) ); } else { //cout << " this says that " << poly_ex << " coeff for degree "<< i <<" is " << cf << endl; } } return result; }; struct match_terms : public map_function { lst l_ex; exmap my_factors; lst exp_stack; match_terms(lst& _l) : l_ex(_l) , my_factors() , exp_stack() { } void clear_factors() { my_factors.clear(); } ex operator()(const ex& e) { //cout << " Term is " << e << endl; lst temp_subs; ex _exp_e = e.expand().normal().expand(); //.to_polynomial(temp_subs).normal(); //cout << " expanding term gives us " << _exp_e << endl; //cout << " where " << temp_subs << endl; lst ocurrences; //map< ex , ex > factors; //factors.clear(); coeff_multi_helper( _exp_e , 0 ); dump_factors(); return e; } void dump_factors() { exmap::iterator it; cout << "actual coefficient decomposition: " << endl; for (it = my_factors.begin() ; it != my_factors.end() ; it++) { cout << " coefficients of term : " << it->first << endl; cout << " ::-->> " << it->second << endl; } cout << "no other factors " << endl; } void coeff_multi_helper( const ex& e , int var ) { lst cf; ex c_; if (var == l_ex.nops() ) //last frame { //cf = degree_list( e , l_ex.op( var ) ); //cout << "CMH: var "<< var << "="<< l_ex.op( var ) << " expression: " << e << endl; //cout << "CMH: obtained degree coeffs " << cf << endl; cout << "entering last frame " << endl; exmap::iterator it; ex prev_coeff; ex key = 1; for (int j=0; j< exp_stack.nops() ; j++) { key *= pow( l_ex.op( j ) , exp_stack.op( j ) ); } cout << " before factor find key " << key << endl; dump_factors(); it = my_factors.find( key ); if (it == my_factors.end()) { my_factors[ key ] = collect_common_factors( e ); cout << " coeff of " << key << " <---- " << my_factors[ key ] << endl; cout << "factor dump again: " << endl; dump_factors(); } else { cout << "what the hell we do in here? KEY: " << key << " VALUE: "<< it->second << endl; cout << " wait.. key is " << key << " while " << it->first << " but " << endl; if (key == it->first) cout << "considered equal! " << endl; else cout << " they look diferent " << endl; prev_coeff = it->second; my_factors[ key ] = collect_common_factors( e + prev_coeff ); cout << " coeff of " << key << " <---- " << my_factors[ key ] << endl; dump_factors(); } } else { cf = degree_list( e , l_ex.op( var ) ); cout << "CMH: var "<< var << "="<< l_ex.op( var ) << " expression: " << e << endl; cout << "CMH: obtained degree coeffs " << cf << endl; for (int i=0; i < cf.nops() ; i++ ) { c_ = cf.op( i ); exp_stack.append( c_.op(0) ); //push stack coeff_multi_helper( c_.op(1) , var + 1 ); exp_stack.remove_last(); //pop stack } } } }; int main() { cout << latex; symbol x("x"), y("y"); symbol q("q") , p("p"); cout << "now we should be using map " << endl; //rord = rord.map( post_tr ); lst l_ex; l_ex = lst( x , y , p ); //l_ex = lst( F(q) , F(q).diff(q) , F(q).diff(q,2) , F(q).diff(q,3) , F(q).diff(q,4) ); match_terms mtch( l_ex ); ex rord = pow( x - y , 3 ) + pow( y - p , 3 ) + 3*p*pow( y , 2 ); rord = mtch( rord.expand() ); }; __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas
On Tue, Jan 30, 2007 at 04:43:33PM -0300, Charlls Quarra wrote:
The other problem i had with this program is that there wasnt a builtin way to get coefficients from multivariate polynomials [snipped]
I use this function to convert a polynomial into an exmap (std::map<ex, ex, ex_is_less>): exmap vargs_collect(const ex& e, const lst& l) { ex x = e.expand(); exmap cmap; cmap[_ex1] = _ex0; for (const_iterator xi=x.begin(); xi!=x.end(); ++xi) { ex key = _ex1; ex pre_coeff = *xi; for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) { int cexp = pre_coeff.degree(*li); pre_coeff = pre_coeff.coeff(*li, cexp); key *= pow(*li, cexp); } exmap::iterator ci = cmap.find(key); if (ci != cmap.end()) ci->second += pre_coeff; else cmap.insert(exmap::value_type(key, pre_coeff)); } return cmap; } After such a conversion getting coefficients is trivial: ex test = pow(x+y+z, 10); exmap tm = vargs_collect(test, lst(x, y)); cout << tm[pow(x,7)*y] << endl; // 360*z^2 Best regards, Alexei -- All science is either physics or stamp collecting.
participants (2)
-
Charlls Quarra
-
varg＠theor.jinr.ru