Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
812 views
in Technique[技术] by (71.8m points)

python - What does r() function mean in the return value of SymPy's dsolve?

I want to evaluate the value of phi(+oo) where phi(xi) is the solution of ODE

Eq(Derivative(phi(xi), (xi, 2)), (-K + xi**2)*phi(xi))

and K is a known real variable. By dsolve, I got the solution:

Eq(phi(xi), -K*xi**5*r(3)/20 + C2*(K**2*xi**4/24 - K*xi**2/2 + xi**4/12 + 1) + C1*xi*(xi**4/20 + 1) + O(xi**6))

with an unknown function r() in the first term on the right-hand side. Here is my code:

import numpy as np
import matplotlib.pyplot as plt
import sympy
from sympy import I, pi, oo

sympy.init_printing()

def apply_ics(sol, ics, x, known_params):
    """
    Apply the initial conditions (ics), given as a dictionary on
    the form ics = {y(0): y0, y(x).diff(x).subs(x, 0): yp0, ...},
    to the solution of the ODE with independent variable x.
    The undetermined integration constants C1, C2, ... are extracted
    from the free symbols of the ODE solution, excluding symbols in
    the known_params list.
    """

    free_params = sol.free_symbols - set(known_params)
    eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics)
            for n in range(len(ics))]
    sol_params = sympy.solve(eqs, free_params)

    return sol.subs(sol_params)

K = sympy.Symbol('K', positive = True)
xi = sympy.Symbol('xi',real = True)
phi = sympy.Function('phi')
ode = sympy.Eq( phi(xi).diff(xi, 2), (xi**2-K)*phi(xi))

ode_sol = sympy.dsolve(ode)
ics = { phi(0):1, phi(xi).diff(xi).subs(xi,0): 0}
phi_xi_sol = apply_ics(ode_sol, ics, xi, [K])

Where ode_sol is the solution, phi_xi_sol is the solution after initial conditions are applied. Since r() is undefined in NumPy I can't evaluate the results by

for g in [0.9, 0.95, 1, 1.05, 1.2]:
    phi_xi = sympy.lambdify(xi, phi_xi_sol.rhs.subs({K:g}), 'numpy')

Does anyone know what this function r() mean and how should I deal with it?

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Answer

0 votes
by (71.8m points)

As visible in the form of the result, the solver falls back to a power series solution (instead of searching the solution in terms of parabolic cylinder functions as WolframAlpha does).

So let's set phi(xi)=sum a[k]*xi^k leading to the coefficient equations (using a[k]=0 for k<0)

(k+2)(k+1)a[k+2] = -K*a[k] + a[k-2]

a[0] = C2  
a[1] = C1  
a[2] = -K/2*C2
a[3] = -K/6*C1
a[4] = (K^2/2 + 1)/12*C2
a[5] = (K^2/6 + 1)/20*C1

Inserting that the the power series solution should have been

C2*(1-K/2*xi**2+(K**2/24+1/12)*xi**4) + C1*xi*(1-K/6*xi**2+(K/120+1/20)*xi**4) + O(xi**6)

Comparing with the sympy solution, all terms containing both C1 and K are missing, especially the missing degree 3 term is not explainable. It seems that the solution process was prematurely ended, or some equation transformation was not correctly reversed.

Please note that the ODE solver routines in sympy are experimental and rudimentary. Also, the power series solution gives only valid information for small values of xi, there is no way to derive any exact value for the limit at +oo.


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...