Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

648

649

650

651

652

653

654

655

656

657

658

659

660

661

662

663

664

665

666

667

668

669

670

671

672

673

674

675

676

677

678

679

680

681

682

683

684

685

686

687

688

689

690

691

692

693

694

695

696

697

698

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719

720

721

722

723

724

725

726

727

728

729

730

731

732

733

734

735

736

737

738

739

740

741

742

743

744

745

746

747

748

749

750

751

752

753

754

755

756

757

758

759

760

761

762

763

764

765

766

767

768

769

770

771

772

773

774

775

776

777

778

779

780

781

782

783

784

785

786

787

788

789

790

791

792

793

794

795

796

797

798

799

800

801

802

803

804

805

806

807

808

809

810

811

812

813

814

815

816

817

818

819

820

821

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

838

839

840

841

842

843

844

845

846

847

848

849

850

851

852

853

854

855

856

857

858

859

860

861

862

863

864

865

866

867

868

869

870

871

872

873

874

875

876

877

878

879

880

881

882

883

884

885

886

887

888

889

890

891

892

893

894

895

896

897

898

899

900

901

902

903

904

905

906

907

908

909

910

911

912

913

914

915

916

917

918

919

920

921

922

923

924

925

926

927

928

929

930

931

932

933

934

935

936

937

938

939

940

941

942

943

944

945

946

947

948

949

950

951

952

953

954

955

956

957

958

959

960

961

962

963

964

965

966

967

968

969

970

971

972

973

974

975

976

977

978

979

980

981

982

983

984

985

986

987

988

989

990

991

992

993

994

995

996

997

998

""" 

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

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

infrastructure remains the same. 

 

**Functions in this module** 

 

    These are the user functions in this module: 

 

    - pdsolve()     - Solves PDE's 

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

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

                       additive or multiplicative separation approach. 

 

    These are the helper functions in this module: 

 

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

    - pde_separate_mul() - Helper function for searching multiplicative 

                           separable solutions. 

 

**Currently implemented solver methods** 

 

The following methods are implemented for solving partial differential 

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

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

 

  - 1st order linear homogeneous partial differential equations 

    with constant coefficients. 

  - 1st order linear general partial differential equations 

    with constant coefficients. 

  - 1st order linear partial differential equations with 

    variable coefficients. 

 

""" 

from __future__ import print_function, division 

 

from sympy.simplify import simplify 

from sympy.core import Add, S 

from sympy.core.compatibility import (reduce, combinations_with_replacement, 

    is_sequence, range) 

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

from sympy.core.relational import Equality, Eq 

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

from sympy.functions import exp 

from sympy.integrals.integrals import Integral 

from sympy.utilities.iterables import has_dups 

 

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

from sympy.solvers.solvers import solve 

from sympy.simplify.simplify import collect 

import operator 

 

allhints = ( 

    "1st_linear_constant_coeff_homogeneous", 

    "1st_linear_constant_coeff", 

    "1st_linear_constant_coeff_Integral", 

    "1st_linear_variable_coeff" 

    ) 

 

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

    """ 

    Solves any (supported) kind of partial differential equation. 

 

    **Usage** 

 

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

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

 

    **Details** 

 

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

            the pde docstring for supported methods).  This can either 

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

            equal to 0. 

 

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

            variable make up the partial differential equation. In many 

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

            (and an error raised if it couldn't be detected). 

 

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

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

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

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

            more options that you can use for hint. 

 

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

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

            to be F. 

 

    **Hints** 

 

        Aside from the various solving methods, there are also some 

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

 

        "default": 

                This uses whatever hint is returned first by 

                classify_pde(). This is the default argument to 

                pdsolve(). 

 

        "all": 

                To make pdsolve apply all relevant classification hints, 

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

                dictionary of hint:solution terms.  If a hint causes 

                pdsolve to raise the NotImplementedError, value of that 

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

                dictionary will also include some special keys: 

 

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

                  deutils.py 

                - default: The solution that would be returned by 

                  default.  This is the one produced by the hint that 

                  appears first in the tuple returned by classify_pde(). 

 

        "all_Integral": 

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

                corresponding "_Integral" hint, it only returns the 

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

                pdsolve() to hang because of a difficult or impossible 

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

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

 

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

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

 

    **Tips** 

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

 

            >>> from sympy import Function, Derivative 

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

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

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

            >>> fx = Derivative(f, x) 

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

            >>> fy = Derivative(f, y) 

 

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

          examples for how to use pdsolve(). 

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

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

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

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

          specific hint 

 

 

    Examples 

    ======== 

 

    >>> from sympy.solvers.pde import pdsolve 

    >>> from sympy import Function, diff, Eq 

    >>> from sympy.abc import x, y 

    >>> f = Function('f') 

    >>> u = f(x, y) 

    >>> ux = u.diff(x) 

    >>> uy = u.diff(y) 

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

    >>> pdsolve(eq) 

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

 

    """ 

 

    given_hint = hint  # hint given by the user. 

 

    if not solvefun: 

        solvefun = Function('F') 

 

    # See the docstring of _desolve for more details. 

    hints = _desolve(eq, func=func, 

        hint=hint, simplify=True, type='pde', **kwargs) 

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

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

 

    if all_: 

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

        # number of hints are added. 

        pdedict = {} 

        failed_hints = {} 

        gethints = classify_pde(eq, dict=True) 

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

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

        for hint in hints: 

            try: 

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

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

            except NotImplementedError as detail: 

                failed_hints[hint] = detail 

            else: 

                pdedict[hint] = rv 

        pdedict.update(failed_hints) 

        return pdedict 

 

    else: 

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

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

 

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

    """Helper function of pdsolve that calls the respective 

    pde functions to solve for the partial differential 

    equations. This minimises the computation in 

    calling _desolve multiple times. 

    """ 

 

    if hint.endswith("_Integral"): 

        solvefunc = globals()[ 

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

    else: 

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

    return _handle_Integral(solvefunc(eq, func, order, 

        match, solvefun), func, order, hint) 

 

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

    r""" 

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

 

    Simplifies the integral mainly using doit() 

    """ 

    if hint.endswith("_Integral"): 

        return expr 

 

    elif hint == "1st_linear_constant_coeff": 

        return simplify(expr.doit()) 

 

    else: 

        return expr 

 

 

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

    """ 

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

 

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

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

    classifications near the beginning of the list will produce 

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

    always exceptions.  To make pdsolve use a different classification, 

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

    docstring for different meta-hints you can use. 

 

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

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

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

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

 

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

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

 

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

    supported hints that can be returned from classify_pde. 

 

 

    Examples 

    ======== 

 

    >>> from sympy.solvers.pde import classify_pde 

    >>> from sympy import Function, diff, Eq 

    >>> from sympy.abc import x, y 

    >>> f = Function('f') 

    >>> u = f(x, y) 

    >>> ux = u.diff(x) 

    >>> uy = u.diff(y) 

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

    >>> classify_pde(eq) 

    ('1st_linear_constant_coeff_homogeneous',) 

    """ 

 

    prep = kwargs.pop('prep', True) 

 

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

        raise NotImplementedError("Right now only partial " 

            "differential equations of two variables are supported") 

 

    if prep or func is None: 

        prep, func_ = _preprocess(eq, func) 

        if func is None: 

            func = func_ 

 

    if isinstance(eq, Equality): 

        if eq.rhs != 0: 

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

        eq = eq.lhs 

 

    f = func.func 

    x = func.args[0] 

    y = func.args[1] 

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

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

 

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

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

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

    # individual variables. 

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

 

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

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

    matching_hints = {'order': order} 

 

    if not order: 

        if dict: 

            matching_hints["default"] = None 

            return matching_hints 

        else: 

            return () 

 

    eq = expand(eq) 

 

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

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

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

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

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

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

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

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

    reduced_eq = None 

    if eq.is_Add: 

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

        dummyvar = var.copy() 

        power = None 

        for i in var: 

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

            if coeff != 1: 

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

                if match and match[a]: 

                    power = match[n] 

                    dummyvar.remove(i) 

                    break 

            dummyvar.remove(i) 

        for i in dummyvar: 

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

            if coeff != 1: 

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

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

                    power = match[n] 

        if power: 

            den = f(x,y)**power 

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

    if not reduced_eq: 

        reduced_eq = eq 

 

    if order == 1: 

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

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

        if r: 

            if not r[e]: 

                ## Linear first-order homogeneous partial-differential 

                ## equation with constant coefficients 

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

                matching_hints["1st_linear_constant_coeff_homogeneous"] = r 

            else: 

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

                    ## Linear first-order general partial-differential 

                    ## equation with constant coefficients 

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

                    matching_hints["1st_linear_constant_coeff"] = r 

                    matching_hints[ 

                        "1st_linear_constant_coeff_Integral"] = r 

 

        else: 

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

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

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

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

            if r: 

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

                matching_hints["1st_linear_variable_coeff"] = r 

 

    # Order keys based on allhints. 

    retlist = [] 

    for i in allhints: 

        if i in matching_hints: 

            retlist.append(i) 

 

    if dict: 

        # Dictionaries are ordered arbitrarily, so make note of which 

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

        matching_hints["default"] = None 

        matching_hints["ordered_hints"] = tuple(retlist) 

        for i in allhints: 

            if i in matching_hints: 

                matching_hints["default"] = i 

                break 

        return matching_hints 

    else: 

        return tuple(retlist) 

 

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

    """ 

    Checks if the given solution satisfies the partial differential 

    equation. 

 

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

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

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

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

    function _preprocess from deutils is used to identify the function. 

 

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

    used to return the result for each solution. 

 

    The following methods are currently being implemented to check if the 

    solution satisfies the PDE: 

 

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

           solution hasn't been solved for f, then it will solve for f 

           provided solve_for_func hasn't been set to False. 

 

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

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

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

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

 

    Examples 

    ======== 

 

    >>> from sympy import Function, symbols, diff 

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

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

    >>> f = Function('f') 

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

    >>> sol = pdsolve(eq) 

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

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

    >>> checkpdesol(eq, sol) 

    (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)) 

    """ 

 

    # Converting the pde into an equation 

    if not isinstance(pde, Equality): 

        pde = Eq(pde, 0) 

 

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

    if func is None: 

        try: 

            _, func = _preprocess(pde.lhs) 

        except ValueError: 

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

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

            funcs = set().union(funcs) 

            if len(funcs) != 1: 

                raise ValueError( 

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

            func = funcs.pop() 

 

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

    # then return a list or set of tuples. 

    if is_sequence(sol, set): 

        return type(sol)([checkpdesol(pde, i, solve_for_func=solve_for_func) for i in sol]) 

 

    # Convert solution into an equation 

    if not isinstance(sol, Equality): 

        sol = Eq(func, sol) 

 

    # Try solving for the function 

    if solve_for_func and not (sol.lhs == func and not sol.rhs.has(func)) and not \ 

            (sol.rhs == func and not sol.lhs.has(func)): 

        try: 

            solved = solve(sol, func) 

            if not solved: 

                raise NotImplementedError 

        except NotImplementedError: 

            pass 

        else: 

            if len(solved) == 1: 

                result = checkpdesol(pde, Eq(func, solved[0]), 

                    order=order, solve_for_func=False) 

            else: 

                result = checkpdesol(pde, [Eq(func, t) for t in solved], 

                order=order, solve_for_func=False) 

 

    # The first method includes direct substitution of the solution in 

    # the PDE and simplifying. 

    pde = pde.lhs - pde.rhs 

    if sol.lhs == func: 

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

    elif sol.rhs == func: 

        s = pde.subs(func, sol.lhs).doit() 

    if s: 

        ss = simplify(s) 

        if ss: 

            return False, ss 

        else: 

            return True, 0 

    else: 

        return True, 0 

 

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

    r""" 

    Solves a first order linear homogeneous 

    partial differential equation with constant coefficients. 

 

    The general form of this partial differential equation is 

 

    .. math:: a \frac{df(x,y)}{dx} + b \frac{df(x,y)}{dy} + c f(x,y) = 0 

 

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

 

    The general solution is of the form:: 

 

        >>> from sympy.solvers import pdsolve 

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

        >>> from sympy import Function, pprint 

        >>> f = Function('f') 

        >>> u = f(x,y) 

        >>> ux = u.diff(x) 

        >>> uy = u.diff(y) 

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

        >>> pprint(genform) 

          d               d 

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

          dx              dy 

 

        >>> pprint(pdsolve(genform)) 

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

                                 --------------- 

                                      2    2 

                                     a  + b 

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

 

    Examples 

    ======== 

 

    >>> from sympy.solvers.pde import ( 

    ... pde_1st_linear_constant_coeff_homogeneous) 

    >>> from sympy import pdsolve 

    >>> from sympy import Function, diff, pprint 

    >>> from sympy.abc import x,y 

    >>> f = Function('f') 

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

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

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

                          x   y 

                        - - - - 

                          2   2 

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

 

    References 

    ========== 

 

    - Viktor Grigoryan, "Partial Differential Equations" 

      Math 124A - Fall 2010, pp.7 

 

    """ 

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

    # two variables are implemented. Once there is support for 

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

 

    f = func.func 

    x = func.args[0] 

    y = func.args[1] 

    b = match[match['b']] 

    c = match[match['c']] 

    d = match[match['d']] 

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

 

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

    r""" 

    Solves a first order linear partial differential equation 

    with constant coefficients. 

 

    The general form of this partial differential equation is 

 

    .. math:: a \frac{df(x,y)}{dx} + b \frac{df(x,y)}{dy} + c f(x,y) = G(x,y) 

 

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

    function in `x` and `y`. 

 

    The general solution of the PDE is:: 

 

        >>> from sympy.solvers import pdsolve 

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

        >>> from sympy import Function, pprint 

        >>> f = Function('f') 

        >>> G = Function('G') 

        >>> u = f(x,y) 

        >>> ux = u.diff(x) 

        >>> uy = u.diff(y) 

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

        >>> pprint(genform) 

                  d               d 

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

                  dx              dy 

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

                  //          b*x + c*y                                             \ 

                  ||              /                                                 | 

                  ||             |                                                  | 

                  ||             |                                       a*xi       | 

                  ||             |                                     -------      | 

                  ||             |                                      2    2      | 

                  ||             |      /b*xi + c*eta  -b*eta + c*xi\  b  + c       | 

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

                  ||             |      |   2    2         2    2   |               | 

                  ||             |      \  b  + c         b  + c    /               | 

                  ||             |                                                  | 

                  ||            /                                                   | 

                  ||                                                                | 

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

                  ||                                  2    2                        | 

                  \\                                 b  + c                         / 

        <BLANKLINE> 

                \| 

                || 

                || 

                || 

                || 

                || 

                || 

                || 

                || 

          -a*xi || 

         -------|| 

          2    2|| 

         b  + c || 

        e       || 

                || 

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

 

 

    Examples 

    ======== 

 

    >>> from sympy.solvers.pde import pdsolve 

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

    >>> from sympy.abc import x,y 

    >>> f = Function('f') 

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

    >>> pdsolve(eq) 

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

 

    References 

    ========== 

 

    - Viktor Grigoryan, "Partial Differential Equations" 

      Math 124A - Fall 2010, pp.7 

 

    """ 

 

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

    # two variables are implemented. Once there is support for 

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

    xi, eta = symbols("xi eta") 

    f = func.func 

    x = func.args[0] 

    y = func.args[1] 

    b = match[match['b']] 

    c = match[match['c']] 

    d = match[match['d']] 

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

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

    functerm = solvefun(eta) 

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

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

    # doit() should be done in _handle_Integral. 

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

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

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

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

 

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

    r""" 

    Solves a first order linear partial differential equation 

    with variable coefficients. The general form of this partial differential equation is 

 

    .. math:: a(x, y) \frac{df(x, y)}{dx} + a(x, y) \frac{df(x, y)}{dy} 

                + c(x, y) f(x, y) - G(x, y) 

 

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

    in `x` and `y`. This PDE is converted into an ODE by making the following transformation. 

 

    1] `\xi` as `x` 

 

    2] `\eta` as the constant in the solution to the differential equation 

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

 

    Making the following substitutions reduces it to the linear ODE 

 

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

 

    which can be solved using dsolve. 

 

    The general form of this PDE is:: 

 

        >>> from sympy.solvers.pde import pdsolve 

        >>> from sympy.abc import x, y 

        >>> from sympy import Function, pprint 

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

        >>> u = f(x,y) 

        >>> ux = u.diff(x) 

        >>> uy = u.diff(y) 

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

        >>> pprint(genform) 

                                             d                     d 

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

                                             dx                    dy 

 

    Examples 

    ======== 

 

    >>> from sympy.solvers.pde import pdsolve 

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

    >>> from sympy.abc import x,y 

    >>> f = Function('f') 

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

    >>> pdsolve(eq) 

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

 

    References 

    ========== 

 

    - Viktor Grigoryan, "Partial Differential Equations" 

      Math 124A - Fall 2010, pp.7 

 

    """ 

    from sympy.integrals.integrals import integrate 

    from sympy.solvers.ode import dsolve 

 

    xi, eta = symbols("xi eta") 

    f = func.func 

    x = func.args[0] 

    y = func.args[1] 

    b = match[match['b']] 

    c = match[match['c']] 

    d = match[match['d']] 

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

 

 

    if not d: 

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

         if not (b and c): 

            if c: 

                try: 

                    tsol = integrate(e/c, y) 

                except NotImplementedError: 

                    raise NotImplementedError("Unable to find a solution" 

                        " due to inability of integrate") 

                else: 

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

            if b: 

                try: 

                    tsol = integrate(e/b, x) 

                except NotImplementedError: 

                    raise NotImplementedError("Unable to find a solution" 

                        " due to inability of integrate") 

                else: 

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

 

    if not c: 

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

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

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

        sol = dsolve(plode, f(x)) 

        syms = sol.free_symbols - plode.free_symbols - set([x, y]) 

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

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

 

    if not b: 

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

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

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

        sol = dsolve(plode, f(y)) 

        syms = sol.free_symbols - plode.free_symbols - set([x, y]) 

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

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

 

    dummy = Function('d') 

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

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

    if isinstance(sol, list): 

        sol = sol[0] 

    solsym = sol.free_symbols - h.free_symbols - set([x, y]) 

    if len(solsym) == 1: 

        solsym = solsym.pop() 

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

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

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

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

        if isinstance(final, list): 

            final = final[0] 

        finsyms = final.free_symbols - deq.free_symbols - set([x, y]) 

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

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

 

    else: 

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

            " to inability of constantsimp") 

 

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

    r""" 

    Helper function to replace constants by functions in 1st_linear_variable_coeff 

    """ 

    eta = Symbol("eta") 

    if len(syms) == 1: 

        sym = syms.pop() 

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

 

    else: 

        fname = func.__name__ 

        for key, sym in enumerate(syms): 

            tempfun = Function(fname + str(key)) 

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

 

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

 

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

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

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

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

    equation than the others. 

 

    :param eq: Partial differential equation 

 

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

 

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

 

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

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

        default. 

 

    Examples 

    ======== 

 

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

    >>> from sympy.abc import x, t 

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

 

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

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

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

 

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

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

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

 

    See Also 

    ======== 

    pde_separate_add, pde_separate_mul 

    """ 

 

    do_add = False 

    if strategy == 'add': 

        do_add = True 

    elif strategy == 'mul': 

        do_add = False 

    else: 

        assert ValueError('Unknown strategy: %s' % strategy) 

 

    if isinstance(eq, Equality): 

        if eq.rhs != 0: 

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

    if eq.rhs != 0: 

        raise ValueError("Value should be 0") 

 

    # Handle arguments 

    orig_args = list(fun.args) 

    subs_args = [] 

    for s in sep: 

        for j in range(0, len(s.args)): 

            subs_args.append(s.args[j]) 

 

    if do_add: 

        functions = reduce(operator.add, sep) 

    else: 

        functions = reduce(operator.mul, sep) 

 

    # Check whether variables match 

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

        raise ValueError("Variable counts do not match") 

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

    if has_dups(subs_args): 

        raise ValueError("Duplicate substitution arguments detected") 

    # Check whether the variables match 

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

        raise ValueError("Arguments do not match") 

 

    # Substitute original function with separated... 

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

 

    # Divide by terms when doing multiplicative separation 

    if not do_add: 

        eq = 0 

        for i in result.args: 

            eq += i/functions 

        result = eq 

 

    svar = subs_args[0] 

    dvar = subs_args[1:] 

    return _separate(result, svar, dvar) 

 

 

def pde_separate_add(eq, fun, sep): 

    """ 

    Helper function for searching additive separable solutions. 

 

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

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

    arguments: 

 

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

 

    Examples 

    ======== 

 

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

    >>> from sympy.abc import x, t 

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

 

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

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

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

 

    """ 

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

 

 

def pde_separate_mul(eq, fun, sep): 

    """ 

    Helper function for searching multiplicative separable solutions. 

 

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

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

    arguments: 

 

    `w(x, y, z) = X(x)*u(y, z)` 

 

    Examples 

    ======== 

 

    >>> from sympy import Function, Eq, pde_separate_mul, Derivative as D 

    >>> from sympy.abc import x, y 

    >>> u, X, Y = map(Function, 'uXY') 

 

    >>> eq = Eq(D(u(x, y), x, 2), D(u(x, y), y, 2)) 

    >>> pde_separate_mul(eq, u(x, y), [X(x), Y(y)]) 

    [Derivative(X(x), x, x)/X(x), Derivative(Y(y), y, y)/Y(y)] 

 

    """ 

    return pde_separate(eq, fun, sep, strategy='mul') 

 

 

def _separate(eq, dep, others): 

    """Separate expression into two parts based on dependencies of variables.""" 

 

    # FIRST PASS 

    # Extract derivatives depending our separable variable... 

    terms = set() 

    for term in eq.args: 

        if term.is_Mul: 

            for i in term.args: 

                if i.is_Derivative and not i.has(*others): 

                    terms.add(term) 

                    continue 

        elif term.is_Derivative and not term.has(*others): 

            terms.add(term) 

    # Find the factor that we need to divide by 

    div = set() 

    for term in terms: 

        ext, sep = term.expand().as_independent(dep) 

        # Failed? 

        if sep.has(*others): 

            return None 

        div.add(ext) 

    # FIXME: Find lcm() of all the divisors and divide with it, instead of 

    # current hack :( 

    # https://github.com/sympy/sympy/issues/4597 

    if len(div) > 0: 

        final = 0 

        for term in eq.args: 

            eqn = 0 

            for i in div: 

                eqn += term / i 

            final += simplify(eqn) 

        eq = final 

 

    # SECOND PASS - separate the derivatives 

    div = set() 

    lhs = rhs = 0 

    for term in eq.args: 

        # Check, whether we have already term with independent variable... 

        if not term.has(*others): 

            lhs += term 

            continue 

        # ...otherwise, try to separate 

        temp, sep = term.expand().as_independent(dep) 

        # Failed? 

        if sep.has(*others): 

            return None 

        # Extract the divisors 

        div.add(sep) 

        rhs -= term.expand() 

    # Do the division 

    fulldiv = reduce(operator.add, div) 

    lhs = simplify(lhs/fulldiv).expand() 

    rhs = simplify(rhs/fulldiv).expand() 

    # ...and check whether we were successful :) 

    if lhs.has(*others) or rhs.has(dep): 

        return None 

    return [lhs, rhs]