Coverage for /usr/lib/python3/dist-packages/sympy/solvers/pde.py: 9%

1""" 

2This module contains pdsolve() and different helper functions that it 

3uses. It is heavily inspired by the ode module and hence the basic 

4infrastructure remains the same. 

5 

6**Functions in this module** 

7 

8 These are the user functions in this module: 

9 

10 - pdsolve() - Solves PDE's 

11 - classify_pde() - Classifies PDEs into possible hints for dsolve(). 

12 - pde_separate() - Separate variables in partial differential equation either by 

13 additive or multiplicative separation approach. 

14 

15 These are the helper functions in this module: 

16 

17 - pde_separate_add() - Helper function for searching additive separable solutions. 

18 - pde_separate_mul() - Helper function for searching multiplicative 

19 separable solutions. 

20 

21**Currently implemented solver methods** 

22 

23The following methods are implemented for solving partial differential 

24equations. See the docstrings of the various pde_hint() functions for 

25more information on each (run help(pde)): 

26 

27 - 1st order linear homogeneous partial differential equations 

28 with constant coefficients. 

29 - 1st order linear general partial differential equations 

30 with constant coefficients. 

31 - 1st order linear partial differential equations with 

32 variable coefficients. 

33 

34""" 

35from functools import reduce 

36 

37from itertools import combinations_with_replacement 

38from sympy.simplify import simplify # type: ignore 

39from sympy.core import Add, S 

40from sympy.core.function import Function, expand, AppliedUndef, Subs 

41from sympy.core.relational import Equality, Eq 

42from sympy.core.symbol import Symbol, Wild, symbols 

43from sympy.functions import exp 

44from sympy.integrals.integrals import Integral, integrate 

45from sympy.utilities.iterables import has_dups, is_sequence 

46from sympy.utilities.misc import filldedent 

47 

48from sympy.solvers.deutils import _preprocess, ode_order, _desolve 

49from sympy.solvers.solvers import solve 

50from sympy.simplify.radsimp import collect 

51 

52import operator 

53 

54 

55allhints = ( 

56 "1st_linear_constant_coeff_homogeneous", 

57 "1st_linear_constant_coeff", 

58 "1st_linear_constant_coeff_Integral", 

59 "1st_linear_variable_coeff" 

60 ) 

61 

62 

63def pdsolve(eq, func=None, hint='default', dict=False, solvefun=None, **kwargs): 

64 """ 

65 Solves any (supported) kind of partial differential equation. 

66 

67 **Usage** 

68 

69 pdsolve(eq, f(x,y), hint) -> Solve partial differential equation 

70 eq for function f(x,y), using method hint. 

71 

72 **Details** 

73 

74 ``eq`` can be any supported partial differential equation (see 

75 the pde docstring for supported methods). This can either 

76 be an Equality, or an expression, which is assumed to be 

77 equal to 0. 

78 

79 ``f(x,y)`` is a function of two variables whose derivatives in that 

80 variable make up the partial differential equation. In many 

81 cases it is not necessary to provide this; it will be autodetected 

82 (and an error raised if it could not be detected). 

83 

84 ``hint`` is the solving method that you want pdsolve to use. Use 

85 classify_pde(eq, f(x,y)) to get all of the possible hints for 

86 a PDE. The default hint, 'default', will use whatever hint 

87 is returned first by classify_pde(). See Hints below for 

88 more options that you can use for hint. 

89 

90 ``solvefun`` is the convention used for arbitrary functions returned 

91 by the PDE solver. If not set by the user, it is set by default 

92 to be F. 

93 

94 **Hints** 

95 

96 Aside from the various solving methods, there are also some 

97 meta-hints that you can pass to pdsolve(): 

98 

99 "default": 

100 This uses whatever hint is returned first by 

101 classify_pde(). This is the default argument to 

102 pdsolve(). 

103 

104 "all": 

105 To make pdsolve apply all relevant classification hints, 

106 use pdsolve(PDE, func, hint="all"). This will return a 

107 dictionary of hint:solution terms. If a hint causes 

108 pdsolve to raise the NotImplementedError, value of that 

109 hint's key will be the exception object raised. The 

110 dictionary will also include some special keys: 

111 

112 - order: The order of the PDE. See also ode_order() in 

113 deutils.py 

114 - default: The solution that would be returned by 

115 default. This is the one produced by the hint that 

116 appears first in the tuple returned by classify_pde(). 

117 

118 "all_Integral": 

119 This is the same as "all", except if a hint also has a 

120 corresponding "_Integral" hint, it only returns the 

121 "_Integral" hint. This is useful if "all" causes 

122 pdsolve() to hang because of a difficult or impossible 

123 integral. This meta-hint will also be much faster than 

124 "all", because integrate() is an expensive routine. 

125 

126 See also the classify_pde() docstring for more info on hints, 

127 and the pde docstring for a list of all supported hints. 

128 

129 **Tips** 

130 - You can declare the derivative of an unknown function this way: 

131 

132 >>> from sympy import Function, Derivative 

133 >>> from sympy.abc import x, y # x and y are the independent variables 

134 >>> f = Function("f")(x, y) # f is a function of x and y 

135 >>> # fx will be the partial derivative of f with respect to x 

136 >>> fx = Derivative(f, x) 

137 >>> # fy will be the partial derivative of f with respect to y 

138 >>> fy = Derivative(f, y) 

139 

140 - See test_pde.py for many tests, which serves also as a set of 

141 examples for how to use pdsolve(). 

142 - pdsolve always returns an Equality class (except for the case 

143 when the hint is "all" or "all_Integral"). Note that it is not possible 

144 to get an explicit solution for f(x, y) as in the case of ODE's 

145 - Do help(pde.pde_hintname) to get help more information on a 

146 specific hint 

147 

148 

149 Examples 

150 ======== 

151 

152 >>> from sympy.solvers.pde import pdsolve 

153 >>> from sympy import Function, Eq 

154 >>> from sympy.abc import x, y 

155 >>> f = Function('f') 

156 >>> u = f(x, y) 

157 >>> ux = u.diff(x) 

158 >>> uy = u.diff(y) 

159 >>> eq = Eq(1 + (2*(ux/u)) + (3*(uy/u)), 0) 

160 >>> pdsolve(eq) 

161 Eq(f(x, y), F(3*x - 2*y)*exp(-2*x/13 - 3*y/13)) 

162 

163 """ 

164 

165 if not solvefun: 

166 solvefun = Function('F') 

167 

168 # See the docstring of _desolve for more details. 

169 hints = _desolve(eq, func=func, hint=hint, simplify=True, 

170 type='pde', **kwargs) 

171 eq = hints.pop('eq', False) 

172 all_ = hints.pop('all', False) 

173 

174 if all_: 

175 # TODO : 'best' hint should be implemented when adequate 

176 # number of hints are added. 

177 pdedict = {} 

178 failed_hints = {} 

179 gethints = classify_pde(eq, dict=True) 

180 pdedict.update({'order': gethints['order'], 

181 'default': gethints['default']}) 

182 for hint in hints: 

183 try: 

184 rv = _helper_simplify(eq, hint, hints[hint]['func'], 

185 hints[hint]['order'], hints[hint][hint], solvefun) 

186 except NotImplementedError as detail: 

187 failed_hints[hint] = detail 

188 else: 

189 pdedict[hint] = rv 

190 pdedict.update(failed_hints) 

191 return pdedict 

192 

193 else: 

194 return _helper_simplify(eq, hints['hint'], hints['func'], 

195 hints['order'], hints[hints['hint']], solvefun) 

196 

197 

198def _helper_simplify(eq, hint, func, order, match, solvefun): 

199 """Helper function of pdsolve that calls the respective 

200 pde functions to solve for the partial differential 

201 equations. This minimizes the computation in 

202 calling _desolve multiple times. 

203 """ 

204 

205 if hint.endswith("_Integral"): 

206 solvefunc = globals()[ 

207 "pde_" + hint[:-len("_Integral")]] 

208 else: 

209 solvefunc = globals()["pde_" + hint] 

210 return _handle_Integral(solvefunc(eq, func, order, 

211 match, solvefun), func, order, hint) 

212 

213 

214def _handle_Integral(expr, func, order, hint): 

215 r""" 

216 Converts a solution with integrals in it into an actual solution. 

217 

218 Simplifies the integral mainly using doit() 

219 """ 

220 if hint.endswith("_Integral"): 

221 return expr 

222 

223 elif hint == "1st_linear_constant_coeff": 

224 return simplify(expr.doit()) 

225 

226 else: 

227 return expr 

228 

229 

230def classify_pde(eq, func=None, dict=False, *, prep=True, **kwargs): 

231 """ 

232 Returns a tuple of possible pdsolve() classifications for a PDE. 

233 

234 The tuple is ordered so that first item is the classification that 

235 pdsolve() uses to solve the PDE by default. In general, 

236 classifications near the beginning of the list will produce 

237 better solutions faster than those near the end, though there are 

238 always exceptions. To make pdsolve use a different classification, 

239 use pdsolve(PDE, func, hint=<classification>). See also the pdsolve() 

240 docstring for different meta-hints you can use. 

241 

242 If ``dict`` is true, classify_pde() will return a dictionary of 

243 hint:match expression terms. This is intended for internal use by 

244 pdsolve(). Note that because dictionaries are ordered arbitrarily, 

245 this will most likely not be in the same order as the tuple. 

246 

247 You can get help on different hints by doing help(pde.pde_hintname), 

248 where hintname is the name of the hint without "_Integral". 

249 

250 See sympy.pde.allhints or the sympy.pde docstring for a list of all 

251 supported hints that can be returned from classify_pde. 

252 

253 

254 Examples 

255 ======== 

256 

257 >>> from sympy.solvers.pde import classify_pde 

258 >>> from sympy import Function, Eq 

259 >>> from sympy.abc import x, y 

260 >>> f = Function('f') 

261 >>> u = f(x, y) 

262 >>> ux = u.diff(x) 

263 >>> uy = u.diff(y) 

264 >>> eq = Eq(1 + (2*(ux/u)) + (3*(uy/u)), 0) 

265 >>> classify_pde(eq) 

266 ('1st_linear_constant_coeff_homogeneous',) 

267 """ 

268 

269 if func and len(func.args) != 2: 

270 raise NotImplementedError("Right now only partial " 

271 "differential equations of two variables are supported") 

272 

273 if prep or func is None: 

274 prep, func_ = _preprocess(eq, func) 

275 if func is None: 

276 func = func_ 

277 

278 if isinstance(eq, Equality): 

279 if eq.rhs != 0: 

280 return classify_pde(eq.lhs - eq.rhs, func) 

281 eq = eq.lhs 

282 

283 f = func.func 

284 x = func.args[0] 

285 y = func.args[1] 

286 fx = f(x,y).diff(x) 

287 fy = f(x,y).diff(y) 

288 

289 # TODO : For now pde.py uses support offered by the ode_order function 

290 # to find the order with respect to a multi-variable function. An 

291 # improvement could be to classify the order of the PDE on the basis of 

292 # individual variables. 

293 order = ode_order(eq, f(x,y)) 

294 

295 # hint:matchdict or hint:(tuple of matchdicts) 

296 # Also will contain "default":<default hint> and "order":order items. 

297 matching_hints = {'order': order} 

298 

299 if not order: 

300 if dict: 

301 matching_hints["default"] = None 

302 return matching_hints 

303 else: 

304 return () 

305 

306 eq = expand(eq) 

307 

308 a = Wild('a', exclude = [f(x,y)]) 

309 b = Wild('b', exclude = [f(x,y), fx, fy, x, y]) 

310 c = Wild('c', exclude = [f(x,y), fx, fy, x, y]) 

311 d = Wild('d', exclude = [f(x,y), fx, fy, x, y]) 

312 e = Wild('e', exclude = [f(x,y), fx, fy]) 

313 n = Wild('n', exclude = [x, y]) 

314 # Try removing the smallest power of f(x,y) 

315 # from the highest partial derivatives of f(x,y) 

316 reduced_eq = None 

317 if eq.is_Add: 

318 var = set(combinations_with_replacement((x,y), order)) 

319 dummyvar = var.copy() 

320 power = None 

321 for i in var: 

322 coeff = eq.coeff(f(x,y).diff(*i)) 

323 if coeff != 1: 

324 match = coeff.match(a*f(x,y)**n) 

325 if match and match[a]: 

326 power = match[n] 

327 dummyvar.remove(i) 

328 break 

329 dummyvar.remove(i) 

330 for i in dummyvar: 

331 coeff = eq.coeff(f(x,y).diff(*i)) 

332 if coeff != 1: 

333 match = coeff.match(a*f(x,y)**n) 

334 if match and match[a] and match[n] < power: 

335 power = match[n] 

336 if power: 

337 den = f(x,y)**power 

338 reduced_eq = Add(*[arg/den for arg in eq.args]) 

339 if not reduced_eq: 

340 reduced_eq = eq 

341 

342 if order == 1: 

343 reduced_eq = collect(reduced_eq, f(x, y)) 

344 r = reduced_eq.match(b*fx + c*fy + d*f(x,y) + e) 

345 if r: 

346 if not r[e]: 

347 ## Linear first-order homogeneous partial-differential 

348 ## equation with constant coefficients 

349 r.update({'b': b, 'c': c, 'd': d}) 

350 matching_hints["1st_linear_constant_coeff_homogeneous"] = r 

351 else: 

352 if r[b]**2 + r[c]**2 != 0: 

353 ## Linear first-order general partial-differential 

354 ## equation with constant coefficients 

355 r.update({'b': b, 'c': c, 'd': d, 'e': e}) 

356 matching_hints["1st_linear_constant_coeff"] = r 

357 matching_hints[ 

358 "1st_linear_constant_coeff_Integral"] = r 

359 

360 else: 

361 b = Wild('b', exclude=[f(x, y), fx, fy]) 

362 c = Wild('c', exclude=[f(x, y), fx, fy]) 

363 d = Wild('d', exclude=[f(x, y), fx, fy]) 

364 r = reduced_eq.match(b*fx + c*fy + d*f(x,y) + e) 

365 if r: 

366 r.update({'b': b, 'c': c, 'd': d, 'e': e}) 

367 matching_hints["1st_linear_variable_coeff"] = r 

368 

369 # Order keys based on allhints. 

370 retlist = [i for i in allhints if i in matching_hints] 

371 

372 if dict: 

373 # Dictionaries are ordered arbitrarily, so make note of which 

374 # hint would come first for pdsolve(). Use an ordered dict in Py 3. 

375 matching_hints["default"] = None 

376 matching_hints["ordered_hints"] = tuple(retlist) 

377 for i in allhints: 

378 if i in matching_hints: 

379 matching_hints["default"] = i 

380 break 

381 return matching_hints 

382 else: 

383 return tuple(retlist) 

384 

385 

386def checkpdesol(pde, sol, func=None, solve_for_func=True): 

387 """ 

388 Checks if the given solution satisfies the partial differential 

389 equation. 

390 

391 pde is the partial differential equation which can be given in the 

392 form of an equation or an expression. sol is the solution for which 

393 the pde is to be checked. This can also be given in an equation or 

394 an expression form. If the function is not provided, the helper 

395 function _preprocess from deutils is used to identify the function. 

396 

397 If a sequence of solutions is passed, the same sort of container will be 

398 used to return the result for each solution. 

399 

400 The following methods are currently being implemented to check if the 

401 solution satisfies the PDE: 

402 

403 1. Directly substitute the solution in the PDE and check. If the 

404 solution has not been solved for f, then it will solve for f 

405 provided solve_for_func has not been set to False. 

406 

407 If the solution satisfies the PDE, then a tuple (True, 0) is returned. 

408 Otherwise a tuple (False, expr) where expr is the value obtained 

409 after substituting the solution in the PDE. However if a known solution 

410 returns False, it may be due to the inability of doit() to simplify it to zero. 

411 

412 Examples 

413 ======== 

414 

415 >>> from sympy import Function, symbols 

416 >>> from sympy.solvers.pde import checkpdesol, pdsolve 

417 >>> x, y = symbols('x y') 

418 >>> f = Function('f') 

419 >>> eq = 2*f(x,y) + 3*f(x,y).diff(x) + 4*f(x,y).diff(y) 

420 >>> sol = pdsolve(eq) 

421 >>> assert checkpdesol(eq, sol)[0] 

422 >>> eq = x*f(x,y) + f(x,y).diff(x) 

423 >>> checkpdesol(eq, sol) 

424 (False, (x*F(4*x - 3*y) - 6*F(4*x - 3*y)/25 + 4*Subs(Derivative(F(_xi_1), _xi_1), _xi_1, 4*x - 3*y))*exp(-6*x/25 - 8*y/25)) 

425 """ 

426 

427 # Converting the pde into an equation 

428 if not isinstance(pde, Equality): 

429 pde = Eq(pde, 0) 

430 

431 # If no function is given, try finding the function present. 

432 if func is None: 

433 try: 

434 _, func = _preprocess(pde.lhs) 

435 except ValueError: 

436 funcs = [s.atoms(AppliedUndef) for s in ( 

437 sol if is_sequence(sol, set) else [sol])] 

438 funcs = set().union(funcs) 

439 if len(funcs) != 1: 

440 raise ValueError( 

441 'must pass func arg to checkpdesol for this case.') 

442 func = funcs.pop() 

443 

444 # If the given solution is in the form of a list or a set 

445 # then return a list or set of tuples. 

446 if is_sequence(sol, set): 

447 return type(sol)([checkpdesol( 

448 pde, i, func=func, 

449 solve_for_func=solve_for_func) for i in sol]) 

450 

451 # Convert solution into an equation 

452 if not isinstance(sol, Equality): 

453 sol = Eq(func, sol) 

454 elif sol.rhs == func: 

455 sol = sol.reversed 

456 

457 # Try solving for the function 

458 solved = sol.lhs == func and not sol.rhs.has(func) 

459 if solve_for_func and not solved: 

460 solved = solve(sol, func) 

461 if solved: 

462 if len(solved) == 1: 

463 return checkpdesol(pde, Eq(func, solved[0]), 

464 func=func, solve_for_func=False) 

465 else: 

466 return checkpdesol(pde, [Eq(func, t) for t in solved], 

467 func=func, solve_for_func=False) 

468 

469 # try direct substitution of the solution into the PDE and simplify 

470 if sol.lhs == func: 

471 pde = pde.lhs - pde.rhs 

472 s = simplify(pde.subs(func, sol.rhs).doit()) 

473 return s is S.Zero, s 

474 

475 raise NotImplementedError(filldedent(''' 

476 Unable to test if %s is a solution to %s.''' % (sol, pde))) 

477 

478 

479 

480def pde_1st_linear_constant_coeff_homogeneous(eq, func, order, match, solvefun): 

481 r""" 

482 Solves a first order linear homogeneous 

483 partial differential equation with constant coefficients. 

484 

485 The general form of this partial differential equation is 

486 

487 .. math:: a \frac{\partial f(x,y)}{\partial x} 

488 + b \frac{\partial f(x,y)}{\partial y} + c f(x,y) = 0 

489 

490 where `a`, `b` and `c` are constants. 

491 

492 The general solution is of the form: 

493 

494 .. math:: 

495 f(x, y) = F(- a y + b x ) e^{- \frac{c (a x + b y)}{a^2 + b^2}} 

496 

497 and can be found in SymPy with ``pdsolve``:: 

498 

499 >>> from sympy.solvers import pdsolve 

500 >>> from sympy.abc import x, y, a, b, c 

501 >>> from sympy import Function, pprint 

502 >>> f = Function('f') 

503 >>> u = f(x,y) 

504 >>> ux = u.diff(x) 

505 >>> uy = u.diff(y) 

506 >>> genform = a*ux + b*uy + c*u 

507 >>> pprint(genform) 

508 d d 

509 a*--(f(x, y)) + b*--(f(x, y)) + c*f(x, y) 

510 dx dy 

511 

512 >>> pprint(pdsolve(genform)) 

513 -c*(a*x + b*y) 

514 --------------- 

515 2 2 

516 a + b 

517 f(x, y) = F(-a*y + b*x)*e 

518 

519 Examples 

520 ======== 

521 

522 >>> from sympy import pdsolve 

523 >>> from sympy import Function, pprint 

524 >>> from sympy.abc import x,y 

525 >>> f = Function('f') 

526 >>> pdsolve(f(x,y) + f(x,y).diff(x) + f(x,y).diff(y)) 

527 Eq(f(x, y), F(x - y)*exp(-x/2 - y/2)) 

528 >>> pprint(pdsolve(f(x,y) + f(x,y).diff(x) + f(x,y).diff(y))) 

529 x y 

530 - - - - 

531 2 2 

532 f(x, y) = F(x - y)*e 

533 

534 References 

535 ========== 

536 

537 - Viktor Grigoryan, "Partial Differential Equations" 

538 Math 124A - Fall 2010, pp.7 

539 

540 """ 

541 # TODO : For now homogeneous first order linear PDE's having 

542 # two variables are implemented. Once there is support for 

543 # solving systems of ODE's, this can be extended to n variables. 

544 

545 f = func.func 

546 x = func.args[0] 

547 y = func.args[1] 

548 b = match[match['b']] 

549 c = match[match['c']] 

550 d = match[match['d']] 

551 return Eq(f(x,y), exp(-S(d)/(b**2 + c**2)*(b*x + c*y))*solvefun(c*x - b*y)) 

552 

553 

554def pde_1st_linear_constant_coeff(eq, func, order, match, solvefun): 

555 r""" 

556 Solves a first order linear partial differential equation 

557 with constant coefficients. 

558 

559 The general form of this partial differential equation is 

560 

561 .. math:: a \frac{\partial f(x,y)}{\partial x} 

562 + b \frac{\partial f(x,y)}{\partial y} 

563 + c f(x,y) = G(x,y) 

564 

565 where `a`, `b` and `c` are constants and `G(x, y)` can be an arbitrary 

566 function in `x` and `y`. 

567 

568 The general solution of the PDE is: 

569 

570 .. math:: 

571 f(x, y) = \left. \left[F(\eta) + \frac{1}{a^2 + b^2} 

572 \int\limits^{a x + b y} G\left(\frac{a \xi + b \eta}{a^2 + b^2}, 

573 \frac{- a \eta + b \xi}{a^2 + b^2} \right) 

574 e^{\frac{c \xi}{a^2 + b^2}}\, d\xi\right] 

575 e^{- \frac{c \xi}{a^2 + b^2}} 

576 \right|_{\substack{\eta=- a y + b x\\ \xi=a x + b y }}\, , 

577 

578 where `F(\eta)` is an arbitrary single-valued function. The solution 

579 can be found in SymPy with ``pdsolve``:: 

580 

581 >>> from sympy.solvers import pdsolve 

582 >>> from sympy.abc import x, y, a, b, c 

583 >>> from sympy import Function, pprint 

584 >>> f = Function('f') 

585 >>> G = Function('G') 

586 >>> u = f(x,y) 

587 >>> ux = u.diff(x) 

588 >>> uy = u.diff(y) 

589 >>> genform = a*ux + b*uy + c*u - G(x,y) 

590 >>> pprint(genform) 

591 d d 

592 a*--(f(x, y)) + b*--(f(x, y)) + c*f(x, y) - G(x, y) 

593 dx dy 

594 >>> pprint(pdsolve(genform, hint='1st_linear_constant_coeff_Integral')) 

595 // a*x + b*y \ 

596 || / | 

597 || | | 

598 || | c*xi | 

599 || | ------- | 

600 || | 2 2 | 

601 || | /a*xi + b*eta -a*eta + b*xi\ a + b | 

602 || | G|------------, -------------|*e d(xi)| 

603 || | | 2 2 2 2 | | 

604 || | \ a + b a + b / | 

605 || | | 

606 || / | 

607 || | 

608 f(x, y) = ||F(eta) + -------------------------------------------------------|* 

609 || 2 2 | 

610 \\ a + b / 

611 <BLANKLINE> 

612 \| 

613 || 

614 || 

615 || 

616 || 

617 || 

618 || 

619 || 

620 || 

621 -c*xi || 

622 -------|| 

623 2 2|| 

624 a + b || 

625 e || 

626 || 

627 /|eta=-a*y + b*x, xi=a*x + b*y 

628 

629 

630 Examples 

631 ======== 

632 

633 >>> from sympy.solvers.pde import pdsolve 

634 >>> from sympy import Function, pprint, exp 

635 >>> from sympy.abc import x,y 

636 >>> f = Function('f') 

637 >>> eq = -2*f(x,y).diff(x) + 4*f(x,y).diff(y) + 5*f(x,y) - exp(x + 3*y) 

638 >>> pdsolve(eq) 

639 Eq(f(x, y), (F(4*x + 2*y)*exp(x/2) + exp(x + 4*y)/15)*exp(-y)) 

640 

641 References 

642 ========== 

643 

644 - Viktor Grigoryan, "Partial Differential Equations" 

645 Math 124A - Fall 2010, pp.7 

646 

647 """ 

648 

649 # TODO : For now homogeneous first order linear PDE's having 

650 # two variables are implemented. Once there is support for 

651 # solving systems of ODE's, this can be extended to n variables. 

652 xi, eta = symbols("xi eta") 

653 f = func.func 

654 x = func.args[0] 

655 y = func.args[1] 

656 b = match[match['b']] 

657 c = match[match['c']] 

658 d = match[match['d']] 

659 e = -match[match['e']] 

660 expterm = exp(-S(d)/(b**2 + c**2)*xi) 

661 functerm = solvefun(eta) 

662 solvedict = solve((b*x + c*y - xi, c*x - b*y - eta), x, y) 

663 # Integral should remain as it is in terms of xi, 

664 # doit() should be done in _handle_Integral. 

665 genterm = (1/S(b**2 + c**2))*Integral( 

666 (1/expterm*e).subs(solvedict), (xi, b*x + c*y)) 

667 return Eq(f(x,y), Subs(expterm*(functerm + genterm), 

668 (eta, xi), (c*x - b*y, b*x + c*y))) 

669 

670 

671def pde_1st_linear_variable_coeff(eq, func, order, match, solvefun): 

672 r""" 

673 Solves a first order linear partial differential equation 

674 with variable coefficients. The general form of this partial 

675 differential equation is 

676 

677 .. math:: a(x, y) \frac{\partial f(x, y)}{\partial x} 

678 + b(x, y) \frac{\partial f(x, y)}{\partial y} 

679 + c(x, y) f(x, y) = G(x, y) 

680 

681 where `a(x, y)`, `b(x, y)`, `c(x, y)` and `G(x, y)` are arbitrary 

682 functions in `x` and `y`. This PDE is converted into an ODE by 

683 making the following transformation: 

684 

685 1. `\xi` as `x` 

686 

687 2. `\eta` as the constant in the solution to the differential 

688 equation `\frac{dy}{dx} = -\frac{b}{a}` 

689 

690 Making the previous substitutions reduces it to the linear ODE 

691 

692 .. math:: a(\xi, \eta)\frac{du}{d\xi} + c(\xi, \eta)u - G(\xi, \eta) = 0 

693 

694 which can be solved using ``dsolve``. 

695 

696 >>> from sympy.abc import x, y 

697 >>> from sympy import Function, pprint 

698 >>> a, b, c, G, f= [Function(i) for i in ['a', 'b', 'c', 'G', 'f']] 

699 >>> u = f(x,y) 

700 >>> ux = u.diff(x) 

701 >>> uy = u.diff(y) 

702 >>> genform = a(x, y)*u + b(x, y)*ux + c(x, y)*uy - G(x,y) 

703 >>> pprint(genform) 

704 d d 

705 -G(x, y) + a(x, y)*f(x, y) + b(x, y)*--(f(x, y)) + c(x, y)*--(f(x, y)) 

706 dx dy 

707 

708 

709 Examples 

710 ======== 

711 

712 >>> from sympy.solvers.pde import pdsolve 

713 >>> from sympy import Function, pprint 

714 >>> from sympy.abc import x,y 

715 >>> f = Function('f') 

716 >>> eq = x*(u.diff(x)) - y*(u.diff(y)) + y**2*u - y**2 

717 >>> pdsolve(eq) 

718 Eq(f(x, y), F(x*y)*exp(y**2/2) + 1) 

719 

720 References 

721 ========== 

722 

723 - Viktor Grigoryan, "Partial Differential Equations" 

724 Math 124A - Fall 2010, pp.7 

725 

726 """ 

727 from sympy.solvers.ode import dsolve 

728 

729 xi, eta = symbols("xi eta") 

730 f = func.func 

731 x = func.args[0] 

732 y = func.args[1] 

733 b = match[match['b']] 

734 c = match[match['c']] 

735 d = match[match['d']] 

736 e = -match[match['e']] 

737 

738 

739 if not d: 

740 # To deal with cases like b*ux = e or c*uy = e 

741 if not (b and c): 

742 if c: 

743 try: 

744 tsol = integrate(e/c, y) 

745 except NotImplementedError: 

746 raise NotImplementedError("Unable to find a solution" 

747 " due to inability of integrate") 

748 else: 

749 return Eq(f(x,y), solvefun(x) + tsol) 

750 if b: 

751 try: 

752 tsol = integrate(e/b, x) 

753 except NotImplementedError: 

754 raise NotImplementedError("Unable to find a solution" 

755 " due to inability of integrate") 

756 else: 

757 return Eq(f(x,y), solvefun(y) + tsol) 

758 

759 if not c: 

760 # To deal with cases when c is 0, a simpler method is used. 

761 # The PDE reduces to b*(u.diff(x)) + d*u = e, which is a linear ODE in x 

762 plode = f(x).diff(x)*b + d*f(x) - e 

763 sol = dsolve(plode, f(x)) 

764 syms = sol.free_symbols - plode.free_symbols - {x, y} 

765 rhs = _simplify_variable_coeff(sol.rhs, syms, solvefun, y) 

766 return Eq(f(x, y), rhs) 

767 

768 if not b: 

769 # To deal with cases when b is 0, a simpler method is used. 

770 # The PDE reduces to c*(u.diff(y)) + d*u = e, which is a linear ODE in y 

771 plode = f(y).diff(y)*c + d*f(y) - e 

772 sol = dsolve(plode, f(y)) 

773 syms = sol.free_symbols - plode.free_symbols - {x, y} 

774 rhs = _simplify_variable_coeff(sol.rhs, syms, solvefun, x) 

775 return Eq(f(x, y), rhs) 

776 

777 dummy = Function('d') 

778 h = (c/b).subs(y, dummy(x)) 

779 sol = dsolve(dummy(x).diff(x) - h, dummy(x)) 

780 if isinstance(sol, list): 

781 sol = sol[0] 

782 solsym = sol.free_symbols - h.free_symbols - {x, y} 

783 if len(solsym) == 1: 

784 solsym = solsym.pop() 

785 etat = (solve(sol, solsym)[0]).subs(dummy(x), y) 

786 ysub = solve(eta - etat, y)[0] 

787 deq = (b*(f(x).diff(x)) + d*f(x) - e).subs(y, ysub) 

788 final = (dsolve(deq, f(x), hint='1st_linear')).rhs 

789 if isinstance(final, list): 

790 final = final[0] 

791 finsyms = final.free_symbols - deq.free_symbols - {x, y} 

792 rhs = _simplify_variable_coeff(final, finsyms, solvefun, etat) 

793 return Eq(f(x, y), rhs) 

794 

795 else: 

796 raise NotImplementedError("Cannot solve the partial differential equation due" 

797 " to inability of constantsimp") 

798 

799 

800def _simplify_variable_coeff(sol, syms, func, funcarg): 

801 r""" 

802 Helper function to replace constants by functions in 1st_linear_variable_coeff 

803 """ 

804 eta = Symbol("eta") 

805 if len(syms) == 1: 

806 sym = syms.pop() 

807 final = sol.subs(sym, func(funcarg)) 

808 

809 else: 

810 for key, sym in enumerate(syms): 

811 final = sol.subs(sym, func(funcarg)) 

812 

813 return simplify(final.subs(eta, funcarg)) 

814 

815 

816def pde_separate(eq, fun, sep, strategy='mul'): 

817 """Separate variables in partial differential equation either by additive 

818 or multiplicative separation approach. It tries to rewrite an equation so 

819 that one of the specified variables occurs on a different side of the 

820 equation than the others. 

821 

822 :param eq: Partial differential equation 

823 

824 :param fun: Original function F(x, y, z) 

825 

826 :param sep: List of separated functions [X(x), u(y, z)] 

827 

828 :param strategy: Separation strategy. You can choose between additive 

829 separation ('add') and multiplicative separation ('mul') which is 

830 default. 

831 

832 Examples 

833 ======== 

834 

835 >>> from sympy import E, Eq, Function, pde_separate, Derivative as D 

836 >>> from sympy.abc import x, t 

837 >>> u, X, T = map(Function, 'uXT') 

838 

839 >>> eq = Eq(D(u(x, t), x), E**(u(x, t))*D(u(x, t), t)) 

840 >>> pde_separate(eq, u(x, t), [X(x), T(t)], strategy='add') 

841 [exp(-X(x))*Derivative(X(x), x), exp(T(t))*Derivative(T(t), t)] 

842 

843 >>> eq = Eq(D(u(x, t), x, 2), D(u(x, t), t, 2)) 

844 >>> pde_separate(eq, u(x, t), [X(x), T(t)], strategy='mul') 

845 [Derivative(X(x), (x, 2))/X(x), Derivative(T(t), (t, 2))/T(t)] 

846 

847 See Also 

848 ======== 

849 pde_separate_add, pde_separate_mul 

850 """ 

851 

852 do_add = False 

853 if strategy == 'add': 

854 do_add = True 

855 elif strategy == 'mul': 

856 do_add = False 

857 else: 

858 raise ValueError('Unknown strategy: %s' % strategy) 

859 

860 if isinstance(eq, Equality): 

861 if eq.rhs != 0: 

862 return pde_separate(Eq(eq.lhs - eq.rhs, 0), fun, sep, strategy) 

863 else: 

864 return pde_separate(Eq(eq, 0), fun, sep, strategy) 

865 

866 if eq.rhs != 0: 

867 raise ValueError("Value should be 0") 

868 

869 # Handle arguments 

870 orig_args = list(fun.args) 

871 subs_args = [arg for s in sep for arg in s.args] 

872 

873 if do_add: 

874 functions = reduce(operator.add, sep) 

875 else: 

876 functions = reduce(operator.mul, sep) 

877 

878 # Check whether variables match 

879 if len(subs_args) != len(orig_args): 

880 raise ValueError("Variable counts do not match") 

881 # Check for duplicate arguments like [X(x), u(x, y)] 

882 if has_dups(subs_args): 

883 raise ValueError("Duplicate substitution arguments detected") 

884 # Check whether the variables match 

885 if set(orig_args) != set(subs_args): 

886 raise ValueError("Arguments do not match") 

887 

888 # Substitute original function with separated... 

889 result = eq.lhs.subs(fun, functions).doit() 

890 

891 # Divide by terms when doing multiplicative separation 

892 if not do_add: 

893 eq = 0 

894 for i in result.args: 

895 eq += i/functions 

896 result = eq 

897 

898 svar = subs_args[0] 

899 dvar = subs_args[1:] 

900 return _separate(result, svar, dvar) 

901 

902 

903def pde_separate_add(eq, fun, sep): 

904 """ 

905 Helper function for searching additive separable solutions. 

906 

907 Consider an equation of two independent variables x, y and a dependent 

908 variable w, we look for the product of two functions depending on different 

909 arguments: 

910 

911 `w(x, y, z) = X(x) + y(y, z)` 

912 

913 Examples 

914 ======== 

915 

916 >>> from sympy import E, Eq, Function, pde_separate_add, Derivative as D 

917 >>> from sympy.abc import x, t 

918 >>> u, X, T = map(Function, 'uXT') 

919 

920 >>> eq = Eq(D(u(x, t), x), E**(u(x, t))*D(u(x, t), t)) 

921 >>> pde_separate_add(eq, u(x, t), [X(x), T(t)]) 

922 [exp(-X(x))*Derivative(X(x), x), exp(T(t))*Derivative(T(t), t)] 

923 

924 """ 

925 return pde_separate(eq, fun, sep, strategy='add') 

926 

927 

928def pde_separate_mul(eq, fun, sep): 

929 """