Coverage for /usr/lib/python3/dist-packages/sympy/solvers/ode/hypergeometric.py: 9%
171 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1r'''
2This module contains the implementation of the 2nd_hypergeometric hint for
3dsolve. This is an incomplete implementation of the algorithm described in [1].
4The algorithm solves 2nd order linear ODEs of the form
6.. math:: y'' + A(x) y' + B(x) y = 0\text{,}
8where `A` and `B` are rational functions. The algorithm should find any
9solution of the form
11.. math:: y = P(x) _pF_q(..; ..;\frac{\alpha x^k + \beta}{\gamma x^k + \delta})\text{,}
13where pFq is any of 2F1, 1F1 or 0F1 and `P` is an "arbitrary function".
14Currently only the 2F1 case is implemented in SymPy but the other cases are
15described in the paper and could be implemented in future (contributions
16welcome!).
18References
19==========
21.. [1] L. Chan, E.S. Cheb-Terrab, Non-Liouvillian solutions for second order
22 linear ODEs, (2004).
23 https://arxiv.org/abs/math-ph/0402063
24'''
26from sympy.core import S, Pow
27from sympy.core.function import expand
28from sympy.core.relational import Eq
29from sympy.core.symbol import Symbol, Wild
30from sympy.functions import exp, sqrt, hyper
31from sympy.integrals import Integral
32from sympy.polys import roots, gcd
33from sympy.polys.polytools import cancel, factor
34from sympy.simplify import collect, simplify, logcombine # type: ignore
35from sympy.simplify.powsimp import powdenest
36from sympy.solvers.ode.ode import get_numbered_constants
39def match_2nd_hypergeometric(eq, func):
40 x = func.args[0]
41 df = func.diff(x)
42 a3 = Wild('a3', exclude=[func, func.diff(x), func.diff(x, 2)])
43 b3 = Wild('b3', exclude=[func, func.diff(x), func.diff(x, 2)])
44 c3 = Wild('c3', exclude=[func, func.diff(x), func.diff(x, 2)])
45 deq = a3*(func.diff(x, 2)) + b3*df + c3*func
46 r = collect(eq,
47 [func.diff(x, 2), func.diff(x), func]).match(deq)
48 if r:
49 if not all(val.is_polynomial() for val in r.values()):
50 n, d = eq.as_numer_denom()
51 eq = expand(n)
52 r = collect(eq, [func.diff(x, 2), func.diff(x), func]).match(deq)
54 if r and r[a3]!=0:
55 A = cancel(r[b3]/r[a3])
56 B = cancel(r[c3]/r[a3])
57 return [A, B]
58 else:
59 return []
62def equivalence_hypergeometric(A, B, func):
63 # This method for finding the equivalence is only for 2F1 type.
64 # We can extend it for 1F1 and 0F1 type also.
65 x = func.args[0]
67 # making given equation in normal form
68 I1 = factor(cancel(A.diff(x)/2 + A**2/4 - B))
70 # computing shifted invariant(J1) of the equation
71 J1 = factor(cancel(x**2*I1 + S(1)/4))
72 num, dem = J1.as_numer_denom()
73 num = powdenest(expand(num))
74 dem = powdenest(expand(dem))
75 # this function will compute the different powers of variable(x) in J1.
76 # then it will help in finding value of k. k is power of x such that we can express
77 # J1 = x**k * J0(x**k) then all the powers in J0 become integers.
78 def _power_counting(num):
79 _pow = {0}
80 for val in num:
81 if val.has(x):
82 if isinstance(val, Pow) and val.as_base_exp()[0] == x:
83 _pow.add(val.as_base_exp()[1])
84 elif val == x:
85 _pow.add(val.as_base_exp()[1])
86 else:
87 _pow.update(_power_counting(val.args))
88 return _pow
90 pow_num = _power_counting((num, ))
91 pow_dem = _power_counting((dem, ))
92 pow_dem.update(pow_num)
94 _pow = pow_dem
95 k = gcd(_pow)
97 # computing I0 of the given equation
98 I0 = powdenest(simplify(factor(((J1/k**2) - S(1)/4)/((x**k)**2))), force=True)
99 I0 = factor(cancel(powdenest(I0.subs(x, x**(S(1)/k)), force=True)))
101 # Before this point I0, J1 might be functions of e.g. sqrt(x) but replacing
102 # x with x**(1/k) should result in I0 being a rational function of x or
103 # otherwise the hypergeometric solver cannot be used. Note that k can be a
104 # non-integer rational such as 2/7.
105 if not I0.is_rational_function(x):
106 return None
108 num, dem = I0.as_numer_denom()
110 max_num_pow = max(_power_counting((num, )))
111 dem_args = dem.args
112 sing_point = []
113 dem_pow = []
114 # calculating singular point of I0.
115 for arg in dem_args:
116 if arg.has(x):
117 if isinstance(arg, Pow):
118 # (x-a)**n
119 dem_pow.append(arg.as_base_exp()[1])
120 sing_point.append(list(roots(arg.as_base_exp()[0], x).keys())[0])
121 else:
122 # (x-a) type
123 dem_pow.append(arg.as_base_exp()[1])
124 sing_point.append(list(roots(arg, x).keys())[0])
126 dem_pow.sort()
127 # checking if equivalence is exists or not.
129 if equivalence(max_num_pow, dem_pow) == "2F1":
130 return {'I0':I0, 'k':k, 'sing_point':sing_point, 'type':"2F1"}
131 else:
132 return None
135def match_2nd_2F1_hypergeometric(I, k, sing_point, func):
136 x = func.args[0]
137 a = Wild("a")
138 b = Wild("b")
139 c = Wild("c")
140 t = Wild("t")
141 s = Wild("s")
142 r = Wild("r")
143 alpha = Wild("alpha")
144 beta = Wild("beta")
145 gamma = Wild("gamma")
146 delta = Wild("delta")
147 # I0 of the standerd 2F1 equation.
148 I0 = ((a-b+1)*(a-b-1)*x**2 + 2*((1-a-b)*c + 2*a*b)*x + c*(c-2))/(4*x**2*(x-1)**2)
149 if sing_point != [0, 1]:
150 # If singular point is [0, 1] then we have standerd equation.
151 eqs = []
152 sing_eqs = [-beta/alpha, -delta/gamma, (delta-beta)/(alpha-gamma)]
153 # making equations for the finding the mobius transformation
154 for i in range(3):
155 if i<len(sing_point):
156 eqs.append(Eq(sing_eqs[i], sing_point[i]))
157 else:
158 eqs.append(Eq(1/sing_eqs[i], 0))
159 # solving above equations for the mobius transformation
160 _beta = -alpha*sing_point[0]
161 _delta = -gamma*sing_point[1]
162 _gamma = alpha
163 if len(sing_point) == 3:
164 _gamma = (_beta + sing_point[2]*alpha)/(sing_point[2] - sing_point[1])
165 mob = (alpha*x + beta)/(gamma*x + delta)
166 mob = mob.subs(beta, _beta)
167 mob = mob.subs(delta, _delta)
168 mob = mob.subs(gamma, _gamma)
169 mob = cancel(mob)
170 t = (beta - delta*x)/(gamma*x - alpha)
171 t = cancel(((t.subs(beta, _beta)).subs(delta, _delta)).subs(gamma, _gamma))
172 else:
173 mob = x
174 t = x
176 # applying mobius transformation in I to make it into I0.
177 I = I.subs(x, t)
178 I = I*(t.diff(x))**2
179 I = factor(I)
180 dict_I = {x**2:0, x:0, 1:0}
181 I0_num, I0_dem = I0.as_numer_denom()
182 # collecting coeff of (x**2, x), of the standerd equation.
183 # substituting (a-b) = s, (a+b) = r
184 dict_I0 = {x**2:s**2 - 1, x:(2*(1-r)*c + (r+s)*(r-s)), 1:c*(c-2)}
185 # collecting coeff of (x**2, x) from I0 of the given equation.
186 dict_I.update(collect(expand(cancel(I*I0_dem)), [x**2, x], evaluate=False))
187 eqs = []
188 # We are comparing the coeff of powers of different x, for finding the values of
189 # parameters of standerd equation.
190 for key in [x**2, x, 1]:
191 eqs.append(Eq(dict_I[key], dict_I0[key]))
193 # We can have many possible roots for the equation.
194 # I am selecting the root on the basis that when we have
195 # standard equation eq = x*(x-1)*f(x).diff(x, 2) + ((a+b+1)*x-c)*f(x).diff(x) + a*b*f(x)
196 # then root should be a, b, c.
198 _c = 1 - factor(sqrt(1+eqs[2].lhs))
199 if not _c.has(Symbol):
200 _c = min(list(roots(eqs[2], c)))
201 _s = factor(sqrt(eqs[0].lhs + 1))
202 _r = _c - factor(sqrt(_c**2 + _s**2 + eqs[1].lhs - 2*_c))
203 _a = (_r + _s)/2
204 _b = (_r - _s)/2
206 rn = {'a':simplify(_a), 'b':simplify(_b), 'c':simplify(_c), 'k':k, 'mobius':mob, 'type':"2F1"}
207 return rn
210def equivalence(max_num_pow, dem_pow):
211 # this function is made for checking the equivalence with 2F1 type of equation.
212 # max_num_pow is the value of maximum power of x in numerator
213 # and dem_pow is list of powers of different factor of form (a*x b).
214 # reference from table 1 in paper - "Non-Liouvillian solutions for second order
215 # linear ODEs" by L. Chan, E.S. Cheb-Terrab.
216 # We can extend it for 1F1 and 0F1 type also.
218 if max_num_pow == 2:
219 if dem_pow in [[2, 2], [2, 2, 2]]:
220 return "2F1"
221 elif max_num_pow == 1:
222 if dem_pow in [[1, 2, 2], [2, 2, 2], [1, 2], [2, 2]]:
223 return "2F1"
224 elif max_num_pow == 0:
225 if dem_pow in [[1, 1, 2], [2, 2], [1, 2, 2], [1, 1], [2], [1, 2], [2, 2]]:
226 return "2F1"
228 return None
231def get_sol_2F1_hypergeometric(eq, func, match_object):
232 x = func.args[0]
233 from sympy.simplify.hyperexpand import hyperexpand
234 from sympy.polys.polytools import factor
235 C0, C1 = get_numbered_constants(eq, num=2)
236 a = match_object['a']
237 b = match_object['b']
238 c = match_object['c']
239 A = match_object['A']
241 sol = None
243 if c.is_integer == False:
244 sol = C0*hyper([a, b], [c], x) + C1*hyper([a-c+1, b-c+1], [2-c], x)*x**(1-c)
245 elif c == 1:
246 y2 = Integral(exp(Integral((-(a+b+1)*x + c)/(x**2-x), x))/(hyperexpand(hyper([a, b], [c], x))**2), x)*hyper([a, b], [c], x)
247 sol = C0*hyper([a, b], [c], x) + C1*y2
248 elif (c-a-b).is_integer == False:
249 sol = C0*hyper([a, b], [1+a+b-c], 1-x) + C1*hyper([c-a, c-b], [1+c-a-b], 1-x)*(1-x)**(c-a-b)
251 if sol:
252 # applying transformation in the solution
253 subs = match_object['mobius']
254 dtdx = simplify(1/(subs.diff(x)))
255 _B = ((a + b + 1)*x - c).subs(x, subs)*dtdx
256 _B = factor(_B + ((x**2 -x).subs(x, subs))*(dtdx.diff(x)*dtdx))
257 _A = factor((x**2 - x).subs(x, subs)*(dtdx**2))
258 e = exp(logcombine(Integral(cancel(_B/(2*_A)), x), force=True))
259 sol = sol.subs(x, match_object['mobius'])
260 sol = sol.subs(x, x**match_object['k'])
261 e = e.subs(x, x**match_object['k'])
263 if not A.is_zero:
264 e1 = Integral(A/2, x)
265 e1 = exp(logcombine(e1, force=True))
266 sol = cancel((e/e1)*x**((-match_object['k']+1)/2))*sol
267 sol = Eq(func, sol)
268 return sol
270 sol = cancel((e)*x**((-match_object['k']+1)/2))*sol
271 sol = Eq(func, sol)
272 return sol