[GiNaC-list] collecting on similar series expansion
Charlls Quarra
charlls_quarra at yahoo.com.ar
Mon Jan 29 23:27:29 CET 2007
--- Sheplyakov Alexei <varg at 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
More information about the GiNaC-list
mailing list