2

Assume that you have a system of coupled PDEs, such as

1st PDE in F(a,b)

2nd PDE in F(a,b)

Following code is able to solve each one of the PDEs separately:

import numpy as np
import sympy as sp

# definition of variables
a, b = sp.symbols('a b')
f = sp.Function('f')
F = f(a, b)
Fda = F.diff(a)
Fdb = F.diff(b)

# definition of PDEs
eq1 = Fda - 2
eq2 = Fda + Fdb + 2

# solution of separated PDEs
sp.pprint(sp.pdsolve(eq1))
sp.pprint(sp.pdsolve(eq2))

Is it possible to solve the system of PDEs? The syntax would be something like sp.pprint(sp.pdsolve([eq1, eq2])). I have tried [eq1, eq2], {eq1, eq2}, np.array([eq1, eq2]) etc. I had a look at help(sp.pdsolve) and help(sp.pde), but have not yet found a solution.

Karlo
  • 1,630
  • 4
  • 33
  • 57

2 Answers2

3

No, the solution of systems of partial differential equations is not implemented. What actually is implemented:

  1. Solving a 1st order linear PDE with constant coefficients: the general form of solution is known and is hardcoded in the solver; the solver returns it, with given coefficients plugged in.

  2. Solving a 1st order linear PDE with variable coefficients, by converting it to an ODE (known as the method of characteristics). A single PDE only.

Aside: I'm skeptical of symbolic solution of PDE in general, and of systems especially. It's just not a thing that happens outside of carefully composed textbook examples. Either a textbook recipe applies (for textbook problems), or there is a hidden structure to be uncovered with human ingenuity (rare), or there is no symbolic solution to be found.

2

Since your system is separable, it can be solved with dsolve. However, dsolve currently doesn't like things like f(a, b), so you need to manually solve the slices. You also need to manually replace the constants with functions:

>>> fa, fb = symbols('fa fb', cls=Function)
>>> eq1 = fb(a).diff(a) - 2
>>> eq2 = fb(a).diff(a) - fa(b).diff(b) + 2
>>> dsolve(eq1, fb(a))
Eq(fb(a), C1 + 2*a)
>>> fbsol = dsolve(eq1, fb(a)).subs(Symbol("C1"), Function("Ca")(b))
>>> fbsol
Eq(fb(a), 2*a + Ca(b))
>>> eq2.subs(*fbsol.args).doit()
Derivative(fa(b), b) + 4
>>> fasol = dsolve(eq2.subs(*fbsol.args).doit(), fa(b)).subs(Symbol("C1"), Function("Cb")(a))
>>> fasol
Eq(fa(b), -4*b + Cb(a))
>>> fbsol
Eq(fb(a), 2*a + Ca(b))

From here it should be clear that Cb(a) = 2*a + C and Ca(b) = -4*b + C, giving the solution f(a, b) = 2*a - 4*b + C, which you can check satisfies the original.

This is definitely something that pdsolve should be able to do automatically, but not much is implemented there yet.

asmeurer
  • 86,894
  • 26
  • 169
  • 240