On Thu, Aug 30, 2012 at 10:56 AM, Aaron Meurer <[email protected]> wrote: > On Thu, Aug 30, 2012 at 1:28 AM, Ondřej Čertík <[email protected]> > wrote: >> Hi, >> >> Let's say that I have just first derivative (but I would like second >> derivatives too later): >> >> >>>>> f(x).diff(x) >> d >> ──(f(x)) >> dx >> >> and I would like to use the substitution y = 2*x**2 to transform the >> derivative. Mathematically, here is how it works: >> >> df/dx = df/dy * dy/dx = df/dy * 4*x >> >> I can do the following: >> >> >>>>> eq = f(x).diff(x) >>>>> eq.subs(x, sqrt(y/2)) >> ⎛d ⎞│ ___ ___ >> ⎜──(f(x))⎟│ ╲╱ 2 ⋅╲╱ y >> ⎝dx ⎠│x=─────────── >> 2 >> >> it's a little ugly because I have to manually invert the substitution >> and also the derivative is still expressed in terms of "x". I would >> like to get: >> >> f'(y) * 4*x >> >> or >> >> f'(y) * 4*sqrt(y/2) >> >> I was wondering whether there is any support for changing a symbolic >> ODE from one variable to another. > > I don't think there is, but obviously it would be useful. Would this > make more sense as a method on Expr or as a function in ode.py? > >> The confusing notion is that if you write f(x) and f(y), it's not >> clear whether it's the same function, just relabeling the independent >> variable, >> or whether those are two different functions, one is f(x), and the >> other is f(y) = f(2*x**2). > > So the real issue is that y should really be y(x), another function.
Yes. Maybe it's better to do it on an actual problem. Here is my script for radial Schroedinger equation for Coulombic potential: https://gist.github.com/3541041 You can see the output there. The idea is simple. You start from this equation: -(-2*E - 2*Z/r + l*(l + 1)/r**2)*P(r) + Derivative(P(r), r, r) (it's always assumed that this is the left hand side, and the right hand side is equal to 0). You first do the substitution lam = sqrt(-2*E), but that's just cosmetic: -(-2*Z/r + l*(l + 1)/r**2 + lam**2)*P(r) + Derivative(P(r), r, r) Now the big part is that you factor out the asymptotic form for r->0 and r->oo and do the substitution: P(r) = r^{l+1}*e^{-lam*r}*F(r) Here is how I do it (see the script): F = Function("F") eq = eq.subs(P(r), r**(l+1) * exp(-lam*r) * F(r)) # forces the derivative to evaluate and then we simplify the result: eq = eq.doit().simplify() We also divide the whole equation by r^{l+1}*e^{-lam*r}. Here is the result: 2*Z*F(r) - 2*l*lam*F(r) + 2*l*Derivative(F(r), r) - 2*lam*r*Derivative(F(r), r) - 2*lam*F(r) + r*Derivative(F(r), r, r) + 2*Derivative(F(r), r) Then I use .coeff() to collect coefficients: Coefficient F(r): 2*Z - 2*l*lam - 2*lam Coefficient F'(r): 2*l - 2*lam*r + 2 Coefficient F''(r): r The idea now is to convert this into the following form, a so called Kummer's equation that has known solutions using special functions (http://en.wikipedia.org/wiki/Confluent_hypergeometric_function): x F''(x) + (b-x) F'(x) - a F(x) = 0 comparing to the coefficient of F'(r) above, you can see, that the substitution: b = 2(l+1) x = 2*lam*r Might do what we want if we are lucky, but now we need to actually do all the derivatives. I did it on a paper and I got a final equation: x F''(x) + (b-x) F'(x) - a F(x) = 0 with: a = l + 1 - Z/lam b = 2(l+1) x = 2*lam*r The coefficient "a" happens to be this way and everything else nicely cancels. Once you have this in this form, one immediately writes down the exact solution using special functions as well as conditions for the solution to be normalizable, which gives a formula for the energy E (in terms of lam). This is standard stuff, but I want to automate this, so that I can use it for equations where I didn't find the solutions in literature. How can I do the last step, from F(r) into F(x), with the relation x = 2*lam*r? There might even be some automated algorithm to do all the steps above and spit out the final coefficients "a", "b" and "x", but for now I am happy to do it explicitly. Once I can do everything, I'll create a few tests for this, possibly also some helper functions in the ODE module. 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.
