THIS IS A DRAFT! It will appear on the main page once finished!

# Solving wave equation with SymPy

Recently I've been trying to play with wave equation for a guitar/piano string. I wanted to try out different boundary conditions, maybe the inhomogenous versions of equation. (TODO link to manual solution) Anyway, I didn't feel like doing that on paper (especially tedious when you're using something other than $1$ for the various string parameters), and wanted to get something closed form and analytical for performance reasons. Sympy was a natural choice for that... but it turned out to be pretty tricky even for the simplest wave equation! The goal/challenge I set was basically to solve as much as possible without invoking assumptions coming from prior knowledge of the equation. That involves restricting the domains, assuming the specific shape of solution etc.

Hopefully this post gives you an idea of some advanced sympy things you can do, since all the computations using Sympy I've seen so far One could argue that realistically you can't get far with symbolic computations anyway. However, I don't feel like this is a good excuse not to try to improve things and free yourself of the cognitive burden of redoing same symbolic manipulations over and over again.

If you read this post and realise I'm doing more manipulations than I should be doing please do let me know!

# TODO it's a shame it's that hard. TODO something about Nielsen?¶

In [1]:
# use that in Jupiter web client
def mdisplay_md(fmt, *args, **kwargs):
display(Markdown(fmt.format(
*[f'${latex(x)}$' for x in args],
**{k: f'${latex(v)}$' for k, v in kwargs.items()})
))

In [2]:
def as_text(thing):
from IPython.core.interactiveshell import InteractiveShell # type: ignore
plain_formatter = InteractiveShell.instance().display_formatter.formatters['text/plain']
pp = plain_formatter(thing)
lines = pp.splitlines()
return lines

width = len(lines[0])
return [' ' * width for _ in range(height - len(lines))] + lines

# terminal and emacs can't display markdown, so we have to use that as a workaround
def mdisplay_plain(fmt, *args, **kwargs):
import re
from itertools import chain
fargs   = [as_text(a) for a in args]
fkwargs = {k: as_text(v) for k, v in kwargs.items()}

height = max(len(x) for x in chain(fargs, fkwargs.values()))

pargs   = [vpad(a, height) for a in fargs]
pkwargs = {k: vpad(v, height) for k, v in fkwargs.items()}

pos = height // 2

lines = []
for h in range(height):
largs   = [a[h] for a in pargs]
lkwargs = {k: v[h] for k, v in pkwargs.items()}
fstr = fmt if h == pos else re.sub('[^{}]', ' ', fmt)
lines.append(fstr.format(*largs, **lkwargs))
print('\n'.join(lines))

mdisplay = mdisplay_plain

In [3]:
%matplotlib inline
# TODO P hmm try that? %matplotlib notebook . that fucks up export though...
from IPython.display import display, Math, Markdown
from sympy import *
init_printing(use_latex='mathjax')
# TODO P wonder if we can convert to latex by default??
# TODO P I guess the answer is here?? https://docs.sympy.org/latest/_modules/sympy/printing/str.html
# TODO P shit. text has to be displayed as markdown...
# TODO P shit. maybe wrap patterns in text, and everythin else is latex??
# TODO render line numbers
# somehow there is no nice function in sympy to display stuff in the same line...

In [ ]:



Ok, let's define the equation in Sympy now.

In [4]:
from sympy import Derivative as D, Function as F
x, t = symbols('x t', real=True)
u = F('u')(x, t)
eq = Eq(D(u, x, 2), D(u, t, 2))
x_l = 0
L = symbols('L', real=True)
x_r = L

mdisplay("PDE: {}", eq) # TODO velocity??
mdisplay("BCs: $u(x_l, t) = 0$, $u(x_r, t) = 0$, $x_l=${}, $x_r=${}", x_l, x_r)

       2              2
∂              ∂
PDE: ───(u(x, t)) = ───(u(x, t))
2              2
∂x             ∂t
BCs: $u(x_l, t) = 0$, $u(x_r, t) = 0$, $x_l=$0, $x_r=$L


For solving partial differential equations, we've got sympy.solvers.pde.pdsolve. However, there seems to be no way to specify initial/boundary conditions for it. Or well, I guess we'd have to deal with them later. Let's give it a try:

In [5]:
# pdsolve(eq)
# TODO shit! sympy stops at first error...
# TODO find a way to ignore the errors


Huh. It can't handle wave equation?? So, what can it do?

In [6]:
display(pde.allhints)

('1st_linear_constant_coeff_homogeneous',
'1st_linear_constant_coeff',
'1st_linear_constant_coeff_Integral',
'1st_linear_variable_coeff')

Sigh. No second order? That's a bit embarrasing. Ok, turns out we can at least aid with the separation of variables. TODO count hints?

In [7]:
X = F('X')(x)
T = F('T')(t)
[sep_X, sep_T] = pde_separate_mul(eq, u, [X, T])
display((sep_X, sep_T))
# TODO at which point are we actually invoking linearity??

$\displaystyle \left( \frac{\frac{d^{2}}{d x^{2}} X{\left(x \right)}}{X{\left(x \right)}}, \ \frac{\frac{d^{2}}{d t^{2}} T{\left(t \right)}}{T{\left(t \right)}}\right)$

As always with separation of variables, we've have to equate these with some constant. I guess I won't surprise you now by denoting it as $-E$.

In [8]:
E = symbols('E', real=True)
eqX = Eq(sep_X, -E)
eqT = Eq(sep_T, -E)
display((eqX, eqT))

$\displaystyle \left( \frac{\frac{d^{2}}{d x^{2}} X{\left(x \right)}}{X{\left(x \right)}} = - E, \ \frac{\frac{d^{2}}{d t^{2}} T{\left(t \right)}}{T{\left(t \right)}} = - E\right)$

Let's start with solving for $X$. sympy.solvers.ode.dsolve does seem to support boundary conditions! Let's try that. First, let's check whether we can aid Sympy at solving this:

In [9]:
classify_ode(eqX)

Out[9]:
('nth_linear_constant_coeff_homogeneous', '2nd_power_series_ordinary')

Cool, so it should be capable of handling that!

In [10]:
bc_l = Eq(X.subs(x, x_l), 0)
bc_r = Eq(X.subs(x, x_r), 0)
mdisplay("BCS: {}, {}", bc_l, bc_r)
Xsol = dsolve(
eqX,
X,
ics={
bc_l.lhs: bc_l.rhs,
bc_r.lhs: bc_l.rhs,
},
hint='nth_linear_constant_coeff_homogeneous',
)
display(Xsol)

BCS: X(0) = 0, X(L) = 0

$\displaystyle X{\left(x \right)} = 0$

Huh? It certainly is a solution, but we all know there's gotta be more. Let's take a look at the general solution then...

In [11]:
Xsol = dsolve(
eqX,
X,
hint='nth_linear_constant_coeff_homogeneous',
)
display(Xsol)

$\displaystyle X{\left(x \right)} = C_{1} e^{- x \sqrt{- E}} + C_{2} e^{x \sqrt{- E}}$

Right, so looks like the extra solutions are somehow lost in the process. Well, let's try to explicitly solve for the initital conditions.

One problem with that is: how do you know what are the integration variables. I mean, we can see they are $C_1$ and $C_2$, but surely we can't rely on symbols('C_1 C_2') since no one guarantees they would be named same in the next version of sympy. Or you might have had C_1 used up by something else.

First thing I did was to declare consts = Xsol.free_symbols.difference({E, x})). After looking at ode module source code, turned out once could make it even more elegant.

In [12]:
constants = list(Xsol.free_symbols - eqX.free_symbols)
system = [
Eq(Xsol.subs(x, x_l).rhs, 0),
Eq(Xsol.subs(x, x_r).rhs, 0),
]
display(system)

linsolve(system, constants)

$\displaystyle \left[ C_{1} + C_{2} = 0, \ C_{1} e^{- L \sqrt{- E}} + C_{2} e^{L \sqrt{- E}} = 0\right]$
Out[12]:
$\displaystyle \left\{\left( 0, \ 0\right)\right\}$

Trivial solution only? Again? We know there are more, actually infinite number of them, depending on $E$!

After a bit of research turned out there is the solveset module, which is capable of properly solving algebraic equations like $e^{q} = 1$ and actually getting inifinite solution sets. Let's try that:

In [13]:
from sympy.solvers.solveset import linsolve as slinsolve
slinsolve(system, constants)

Out[13]:
$\displaystyle \left\{\left( 0, \ 0\right)\right\}$

Nope.. Perhaps it wants us to specify the parameter we solve against explicitly? We can't use linsolve anymore then since the equation is not linear in $E$, but luckily we've got nonlinsolve!

In [14]:
from sympy.solvers.solveset import nonlinsolve as snonlinsolve
nonlinsolve(system, constants + [E])

Out[14]:
$\displaystyle \left\{\left( 0, \ 0, \ E\right), \left( C_{2}, \ - C_{2}, \ 0\right)\right\}$

No luck. Let's try and debug it by removing the complexity.

A simpler version of that system would be, of course, a system of 1 equation with a parameter, such that discrete values of parameter determine the solutions space.

In [15]:
Y = symbols('Y')
system = [
Y * (exp(sqrt(-E)) - 1)
]
display(*system)

$\displaystyle Y \left(e^{\sqrt{- E}} - 1\right)$
In [16]:
snonlinsolve(system, [Y, E])

Out[16]:
$\displaystyle \left\{\left( 0, \ E\right), \left( Y, \ 0\right)\right\}$

Ok, looks like we're on the right track. It's clearly missing the exponential bit, let's try it in separation.

In [17]:
solveset(exp(sqrt(-E)) - 1)
# TODO need to store captured outputs for persistence

Out[17]:
$\displaystyle \left\{E \mid E \in \mathbb{C} \wedge e^{\sqrt{- E}} - 1 = 0 \right\}$

Right. It's definitely confused. Remember, we defined $E$ to be real, which means that $\sqrt{-E}$ would be pure imaginary and that could be exploited!

And googling immediately yields that

Solveset is designed to be independent of the assumptions on the variable being solved for and instead, uses the domain argument to decide the solver to dispatch the equation to



My gut feeling was that it had something to do with the (multivalued) complex square root. And indeed:

In [18]:
Y = symbols('Y')
solveset(exp(Y) - 1)

Out[18]:
$\displaystyle \left\{2 n i \pi\; |\; n \in \mathbb{Z}\right\}$

Let's experiment for a bit how solveset works.

In [19]:
from sympy.solvers.solveset import solveset_complex
solveset_complex(sqrt(1) + Y, Y)

Out[19]:
$\displaystyle \left\{-1\right\}$

Huh. Surely there should be two solutions, corresponding to $\sqrt 1 = \pm 1$? I suppose it's time to read some docs on how solveset works. This gives us a bit of a hint.

The respective solver now tries to invert the equation using the routines invert_real and invert_complex. These routines are based on the concept of mathematical inverse (though not exactly). It reduces the real/complex valued equation f(x)=y to a set of equations: {g(x)=h1(y),g(x)=h2(y),...,g(x)=hn(y)} where g(x) is a simpler function than f(x).



We're clearly interested at invert_complex code . It looks a bit mad, e.g. here it explicitly excludes trigonometric functions. I wonder why? (TODO raise an issue or something?)

TODO sqrt docs? TODO not sure what that was about...

I figured shit got serious, so I cloned sympy repository and a after a bit, tracked down the problem.

assert solveset_complex(4*x*(1 - a * sqrt(x)), x) == \
FiniteSet(S(0), 1 / a ** 2)



That doesn't look good for us.. there should be an extra solution $-\frac{1}{a^2}$. Git blaming gives us this which tells us the test was commited five years ago. So basically this behaviour has been in sympy for ages and protected by unit test.

I would imagine that the explanation is that multivalued functions like sqrt, asin, etc. are dumbed down and treated as their principal branches.

So in a sense, it's good that we this weird thing $\left\{E \mid E \in \mathbb{C} \wedge e^{\sqrt{E}} - 1 = 0 \right\}$ as a result before, since that at least doesn't assume any specific complex branch. I hope one day I'll get to fix that in Sympy, but for now let's try and at least aid it at converting this thing into the desired solutions set.

Apparently what we need to do is to analyse the subexpressions where $E$ appears under the square root, and then do case analysis on its sign, just as you would do when solving manually. The easy thing to do for the simple wave equation would be just replacing $k = \sqrt{-E}$.

In [20]:
k = symbols('k')
subst = Eq(sqrt(-E), k)

system = [
eq.subs(subst.lhs, subst.rhs) for eq in system
]
Xsol = Xsol.subs(subst.lhs, subst.rhs)
display(system)
display(Xsol)

res = snonlinsolve(system, constants + [k])
mdisplay("Solutions: {}", res)

$\displaystyle \left[ Y \left(e^{k} - 1\right)\right]$
$\displaystyle X{\left(x \right)} = C_{1} e^{- k x} + C_{2} e^{k x}$
Solutions: {(C₂, C₁, 0), (C₂, C₁, {2⋅n⋅ⅈ⋅π | n ∊ ℤ})}


Nice! That looks way better! If you look closer though, you can see that the first solution is contained in the second solution. So we'd like to reduce the complexity and discard the first one.

But if you remember, we want to do it with as little explicit steps as we can. How do we do that? We could try and fully expland the solutions in terms of ImageSet, that is $\{ (C_2, C_1, 2 n i \pi) \mid n \in \mathbb{Z} \}$.

In [21]:
def bindFiniteSet(x):
assert isinstance(x, FiniteSet)
return x.args
def bindImageSet(x):
assert isinstance(x, ImageSet)
return x.lamda, x.base_set
(first, second) = bindFiniteSet(res)
(C1sol, C2sol, kSol) = second
lamda, base = bindImageSet(kSol)
mdisplay("$C_1$ = {}; C2 = {}; k = {}", C1sol, C2sol, kSol)

$C_1$ = C₂; C2 = C₁; k = {2⋅n⋅ⅈ⋅π | n ∊ ℤ}

In [22]:
def bindTuple(x):
assert isinstance(x, Tuple)
return x.args

def flatten(x):
args = bindTuple(x)
args = [x if isinstance(x, Set) else FiniteSet(x) for x in args]
return ProductSet(*args)

def collapse(solutions):
ret = EmptySet()
for sol in solutions.args:
ret = ret.union(flatten(sol))
return ret

collapsed = collapse(res)
display(collapsed)

# TODO need to turn to ImageSet?
# TODO and then we can substitute for each N separately...

$\displaystyle \left\{C_{2}\right\} \times \left\{C_{1}\right\} \times \left\{2 n i \pi\; |\; n \in \mathbb{Z}\right\}$

Now, at the very least it would assert if something goes wrong. Let's actually substitute everything now.

In [23]:
Xsol = Xsol.subs({
constants[0]: C1sol,
constants[1]: C2sol,
})
# only one constant
# [const] = set(constants) - Xsol.free_symbols
# display(const)
A = IndexedBase('A')
B = IndexedBase('B')
n = symbols('n', integer=True)
term = Xsol.rhs.subs({
constants[0]: A[n],
constants[1]: B[n],
k: lamda(n),
})
term = collect(term, [A[n], B[n]]) # for nicer grouping
basis = ImageSet(
Lambda(n, term),
base,
)
mdisplay("Solutions basis: {}", basis) # TODO shit, that doesn't need to be formatted like latex

                 ⎧ 2⋅ⅈ⋅π⋅n⋅x         -2⋅ⅈ⋅π⋅n⋅x             ⎫
Solutions basis: ⎨ℯ         ⋅A[n] + ℯ          ⋅B[n] | n ∊ ℤ⎬
⎩                                          ⎭


Ok, we're getting closer now. One concern is that we've got linearly dependent functions in our basis, the ones corresponding to $n$ and $-n$. I guess it's the nonlinsolve to blame. Can we do anything to eliminate extra solutions in the meantime?

In [24]:
def eliminate_dependency(iset, idx, var):
assert iset.base_set == Integers
prop = simplify(basis.lamda(idx) / basis.lamda(-idx))
print(prop.free_symbols)
assert var not in prop.free_symbols
# now that we verified linear dependency, we can shrink the index space
return ImageSet(iset.lamda, Naturals0)

# basis = eliminate_dependency(basis, n, x)
# display(basis)


Ok, now to find out $B_n$, we're gonna need the general solution. TODO continuus bit -- I suppose another time. TODO actually, integrate over other values of k too! They all result in 0, but still! TODO this might be important when we solve something like schrodinger equation, where the spectrum might have continuous part!

In [25]:
# shit, that didn't work. requires lower and upper bounds
# kinda makes sense, we can't just sum over arbitrary sets...
# but we do want to sum over countable sets!
def sum_iset(iset):
bs = iset.base_set
lamda = iset.lamda
[var] = lamda.variables
return Sum(
lamda(var),
(var, bs._inf, bs._sup),
)
Xsol2 = sum_iset(basis)
display(Xsol2)

$\displaystyle \sum_{n=-\infty}^{\infty} \left(e^{2 i \pi n x} {A}_{n} + e^{- 2 i \pi n x} {B}_{n}\right)$
In [26]:
# TODO FIXME exploiting that A = 0..
def get_func(idx):
return simplify(Xsol2.function.subs(n,idx) / B[idx])

# if we not assume n to be nonnegative it will make integration a bit easier for sympy
# from sympy.assumptions.assume import global_assumptions
# from sympy import Q
m = symbols('m', integer=True)# , nonnegative=True)

# TODO ???? takes a lot to compute this

# integrate(get_func(n) * get_func(m), (x, 0, L))
# TODO wtf??? that's messy

#def sum_map(f, iset, idx):
#    return Sum(
#        f(idx, iset.lamda(idx)),
#        (idx, iset.base_set._inf, iset.base_set._sup),
#    )

# B = IndexedBase('B')

# sol = sum_map(lambda idx, term: Xsol.subs({
#    C1: A[idx],
#    C2: B[idx],
#    k : term,
#}), K, idx=n)

# print('General solution:')
# display(sol)

In [27]:
from sympy.solvers.solveset import invert_complex, invert_real, _solve_radical, _solveset, solveset_complex
QQ = symbols('QQ', real=True)
# snonlinsolve([(Y - 1) * (exp(I * sqrt(YY)) - 1)],[Y, YY])
# TODO FIX -- simplify its output??
from sympy.solvers.solveset import linsolve as slinsolve
display(linsolve(system, constants))
# TODO shit is E even the right name for the constant???

$\displaystyle \emptyset$
In [28]:
from sympy.solvers.solveset import linsolve
from sympy.solvers.ode import ode_nth_linear_constant_coeff_homogeneous as ode_solve
from sympy.core.numbers import IntegerConstant
#
k = symbols('k')
E = k ** 2
# E = symbols('E', real=True)
# eqT = Eq(_eqT + E) # ).canonical
# canonical moves all variables to the left and numbers to the right
Tsol = dsolve(eqT).rhs
# shit, it won't eliminate solutions based on k :(
Xsol = dsolve(
eqX,
X,
#ics={
#    X.subs(x, x_l): 0,
#    X.subs(x, x_r): 0,
#},
hint='nth_linear_constant_coeff_homogeneous',
# returns = 'both',

).rhs

print("General solutions for X and T:")
display((Xsol, Tsol))

# how to specify constants for DE solver??
[C1, C2] = list(Xsol.free_symbols.difference({k, x}))

bcs_l = Eq(Xsol.subs(x, x_l))
bcs_r = Eq(Xsol.subs(x, x_r))

print("BC enforces:")
display((bcs_l, bcs_r))

# TODO blog about that, give examples?
# so, solveset stuff can't handle
[c1_c2] = solve(bcs_l, C1)
# solve(bcs_r.subs(C1, c1_c2))
# [{C2:0},{k:0},{k:iπ}]
# TODO what do we do with C2 = 0 solution???
K = solveset(bcs_r.subs(C1, c1_c2), k)

print("Allowed values of k:")
display(K)

# TODO how to eliminate negative k solutions?
# there must be something wrong with our approach if we have to do that manually. right?
# something unnatural. maybe we want to impose some sort of symmetry first??
# ah, basically separate equation in left and right waves first?

# so, can we eliminate extra constants and instead use the whole range of ks ????

# hmm. maybe the redundancy is coming from separation of variables? e.g. negative and positive energies?

# TODO generate sequence of edits from git commits?...

# basis = imageset(lambda kk: Xsol.subs(k, kk), K)
# print("Solutions basis:")
# display(basis)

# here, sympy retuns something called ImageSet. it basically means that the solutions image of Z under x -> x i pi

# right, how to define infinite coefficients? can it just be a function?
# TODO how to sum something over imageset????

# TODO base_set??

# shit, that didn't work. requires lower and upper bounds
# kinda makes sense, we can't just sum over arbitrary sets...
# but we do want to sum over countable sets!

def sum_map(f, iset, idx):
return Sum(
f(idx, iset.lamda(idx)),
(idx, iset.base_set._inf, iset.base_set._sup),
)

A = IndexedBase('A')
B = IndexedBase('B')
n = symbols('n')

sol = sum_map(lambda idx, term: Xsol.subs({
C1: A[idx],
C2: B[idx],
k : term,
}), K, idx=n)

print('General solution:')
display(sol)

# hmm, it's a bit awkward, we've got quite a bit of redundancy in
# in textbooks, that's violated in an ad-hoc manner at the step (TODO refer to line, with a marker?)
# the proper way to deal with this would be rewriting the solution
# and only then enforce a bc!

# TODO do we need to eliminate negative coefficients??

# imageset(lambda u: u**2, allowed_ks)
#Xsol
# print(allowed_ks.domain)
# ok, now we basically restrict ourselves to considering only specific waves

# TODO pick first N coeffs and do initial transform?
# usol = Xsol.rhs * Tsol.rhs
# linsolve(, [C1, C2])
print("rewrite: C[0] = (A_0 + B_0) / 2; for n>0: C[n] = B_n + A_{-n}; C[-n] = B_{-n} + A_n")
# C = IndexedBase('C')
C = Function('C')
# sol = Sum(C[n] * exp(I * pi * n * x), (n, -oo, oo))

# TODO here, extract the exponent from split solution..
sol = Integral(C(k) * exp(I * k * x), (k, -oo, oo))
print("General sol")
display(sol)
# TODO how to check they are equivalent??
bc_l = sol.subs(x, x_l)
bc_r = sol.subs(x, x_r)
print("BCS:")
display((Eq(bc_l), Eq(bc_r)))
# TODO is there a way to solve that automatically???

def mComm(expr, cls):
assert isinstance(expr, cls)
return expr.args

# TODO is there anything standard in sympy?

def mMul(expr):
return mComm(expr, Mul)

def mIntegral(expr):
assert isinstance(expr, Integral), type(i)
return expr.limits, expr.function

def pullIntegral(expr):
if isinstance(expr, Integral):
return expr
if isinstance(expr, Mul):
# TODO might work for a constant...
parts = [pullIntegral(p) for p in mMul(expr)]
# TODO only one should be integral, rest constants?
consts = []
rest = []
for p in parts:
if isinstance(p, IntegerConstant):
consts.append(p)
else:
rest.append(p)
rest = [pullIntegral(r) for r in rest]
res = None
if len(rest) == 1:
lims, fun = mIntegral(rest[0])
res = Integral(Mul(*consts) * fun, lims)
return res
else:
res = Mul(*consts)
if len(rest) > 0:
res = res * Mul(*rest)
return res
return expr

def asIntegral(expr):
res = pullIntegral(expr)
assert isinstance(res, Integral), type(res)
return res
# if is

def combine_isum(isum):
parts = [asIntegral(p) for p in parts]
# TODO shit, is there something for pulling stuff under integral?
# TODO expression with limits
lims = [i.limits for i in parts]
exprs = [i.function for i in parts]
[lim] = list(set(lims))

# TODO generify for expr with limits

def int2sum(thing):
lims, func = mIntegral(thing)
[q1, _] = mMul(func)
KK = solveset(q1, k)
return KK

# TODO how to bind it to list of integrals? some sort of generic shape??
# TODO right! so if we got independent things, we should test BOTH linear combinations
bcs1 = rcollect(combine_isum(bc_l + bc_r), C(k))
display(bcs1)

bcs2 = rcollect(combine_isum(bc_l - bc_r), C(k))
display(bcs2)

#lims1, func1 = mIntegral(bcs1)
#[qq1, _] = mMul(func1)
#lims2, func2 = mIntegral(bcs2)
#[qq2, _] = mMul(func2)

#from sympy.solvers.solveset import linsolve as slinsolve, nonlinsolve as snonlinsolve
# KK = solveset([qq1, qq2], k)

# KK = slinsolve([qq1, qq2], k)
# display(qq1)
# display(qq2)
# nope.. they were independent and should be solved independently!!

def the(objs):
assert len(objs) > 0, "expected non-empty list"
assert all([o == objs[0] for o in objs]), "not all objects are equal"
return objs[0]

def asLambda1(expr):
assert isinstance(expr, Lambda)
return expr.variables, expr.expr

# TODO https://docs.sympy.org/latest/tutorial/manipulation.html is not too helpful though...
def combine_imagesets(s1, s2):
display(s1)
display(s2)
sets = [s1, s2]
# TODO assert all are image sets? or pull them out
dom = the([s.base_set for s in sets])
# TODO ok, need to  analyse lamdas...
# TODO should be same variable???
v1, e1 = asLambda1(s1.lamda)
v2, e2 = asLambda1(s2.lamda)
# print(cse([e1, e2]))
# assert v1 == v2
# basically, we need to figure out if e1 == e2 for some n? then we can discard one of them!!
# and if only one elem left, we don't need the finiteset?
# ah shit, we need to replace the function!
# display(solve(e1 - e2))
print("GCD")
display(factor_list(e1))
display(factor_list(e2))
display(gcd_terms([e1, e2]))
return imageset(lambda q: FiniteSet(*[s.lamda(q) for s in sets]), dom)

# TODO shit, how to simplify them???
# combine_imagesets()
# fucking hell. why is that so tricky????
# anyway:
display(int2sum(bcs) + int2sum(bcs2))
print("which is same as:")
# TODO again, need some sort of quickcheck for expressions, sets etc
n = symbols('n')
Ns = imageset(Lambda(n, pi * n), Integers)
display(Ns)
# TODO now we actually need int2sum
# something that pick subintegral expression

def change_limits(expr, ns):
# expr.limits
ef = expr.function
[ev] = expr.variables
lamda = ns.lamda
lf = lamda.expr
[lv] = lamda.variables
bs = ns.base_set
return Sum(ef.subs(ev, lamda(lv)), (lv, bs._inf, bs._sup))
# ugh. so dirty..
res2 = change_limits(sol, Ns)
print("After substituting restrictions on wavevectors")
display(res2)

from sympy.integrals.transforms import _fourier_transform

u = symbols('u')
integrate(res2 * exp(-I * pi * u * x), (x, x_l, x_r))

# shit. I guess we run at the fact that expr(-I pi u x) are not basis functions for the problem once again...
# still, shouldn't prevent us from solving. right???

# fucking hell.. why does that work???

# TODO next: initial condition to find Cs?
# basically, here we're going to restrict ourselves? or we can get away with fully symbolic for a while?...
# TODO ok, just multiply all by exp(- i pi k x)
# TODO how to eliminate C_{k+1} vs C_{k-1} ???
# wonder if division by two in classic particle in a box coming from the fact that they have to normalise by 2
# since they are using degenerate solution space

# TODO do final chage of variables from C(pi * n) to some indexed function?..

# TODO shit... but coeffs are not independent.. because originally we had two equations??
# mm. ok, individual things don't have to satisfy bcs!

# when k = 2 pi n, it works for any C
# when it's not, thet that can only work when all of C are 0, due to properties of L^2 space
# TODO could implement solving integral equations in python?...
# basically, C(k) is non-zero on k = 2 pi n and 0 otherwise

# TODO hmm interesting, looks kinda like a discrepancy between pi *n and 2 n pi + pi
# but, actually, I guess this is kinda disguised in the linear combination thing! this form explicitly forbids cosines!

General solutions for X and T:

$\displaystyle \left( C_{1} e^{- x \sqrt{- E}} + C_{2} e^{x \sqrt{- E}}, \ C_{1} e^{- t \sqrt{- E}} + C_{2} e^{t \sqrt{- E}}\right)$
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-28-980207416d34> in <module>
26
27 # how to specify constants for DE solver??
---> 28 [C1, C2] = list(Xsol.free_symbols.difference({k, x}))
29
30 bcs_l = Eq(Xsol.subs(x, x_l))

ValueError: too many values to unpack (expected 2)
In [29]:
# shit. it can't handle restrictions on k :(
from sympy.solvers.solveset import linsolve as slinsolve

A, B, k = symbols('A B k')

slinsolve([
A * exp(sqrt(-k)) + exp(sqrt(k)),
], [A])

# right, it just can't figure out multiple solutions??
# solveset(exp(k) - 1)

# shit, it can't handle eliminating certain k's :(
slinsolve([
A,
0,
#exp(k) - 1,
], [A, B])

E = symbols('E', real=True)
k = symbols('k')

# TODO shit, it won't handle sqrt??
solveset(exp(k) - exp(-k))

# solveset(exp(k) - 1)
solveset(exp(sqrt(k**2)) - 1)

slinsolve([
A + B,
A * exp(- k) + B * exp(k),
],
[A, B]
)

#solve([
#    exp(k) + 1,
#], returns='solveset')

Out[29]:
$\displaystyle \left\{\left( 0, \ 0\right)\right\}$
In [30]:
from sympy.core.rules import Transform
def split(integ):
F = Function('F')
G = Function('G')
H = Function('H')
A, B = symbols('A B')
z = symbols('z')
v = symbols('v')
it = Integral(F(z) * z + -G(z) * z, (z, 0, oo))
display(it)
thing = it.xreplace(Transform(split, lambda i: isinstance(i, Integral)))

left,
right.transform(-z, z),
)
display(rrr)
rrr.subs_expr({
F(z): H(z),
G(-z): H(z),
})
# TODO shit...
# TODO replace function with a partial??

# TODO how to assert that expressions are equivalent???

$\displaystyle \int\limits_{0}^{\infty} \left(z F{\left(z \right)} - z G{\left(z \right)}\right)\, dz$
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-30-a6dca8e88d67> in <module>
12 thing = it.xreplace(Transform(split, lambda i: isinstance(i, Integral)))
13
---> 14 [left, right] = mAdd(thing)
16     left,

NameError: name 'mAdd' is not defined
In [31]:
q, w = symbols('q w')
ee = Eq(q, w)
ee.subs(q, 0)
# display(ee)

Out[31]:
$\displaystyle 0 = w$
In [32]:
# TODO NOTE: the bc equations:
# solve [C1 + C2 =0, C1 * exp(i * k) + C2 * exp(-i * k) =0]
# sage couldn't hande it either
# wolfram alpha did though!
`