1

I recently asked this question about solving a system of nonlinear equations, over on Math SE. I am now attempting to implement the accepted answer using SymPy in Python. I'll begin by briefly recreating the solution here.

Problem.

Suppose I have three receivers, with coordinates given by vectors p1, p2, and p3. Further, I have a transmitter, emitting radio pulses at unknown times t1, t2, t3, t4, ... with known speed c, and which has unknown initial coordinates at time t = t1 given by the vector x0 = <x,y,z> and is moving with an unknown constant velocity given by the vector v. My data consists of the pulse arrival time at each station, denoted air (for arrival of time of signal i at receiver r).

My goal is to determine this velocity. The solution given in the attached answer is as follows:

The position of the transmitter at time t will be given by xo = v(t - t1). Hence, the distance between the transmitter and the rth receiver at the time the ith pulse is emitted will be given by:

||(ti - t1)v + x0 - pr||

Therefore, the time of arrival at the rth receiver is given by:

ti + (1/c) * ||(ti - t1)v + x0 - pr||

Which we know to be air. Therefore:

c**2(air - ti)**2 = ||v||**2(ti - t1)**2 + 2(ti -t1)<v, c0 - pr>

(where <a, b> denotes the scalar/dot product of vectors a and b).

This results in a system of # receivers x # pulses equations.

(see the attached answer for far nicer formatting)

Implementation.

I'm most familiar with the SymPy package, so I'm attempting to implement a solution using its solvers. I'm a bit stuck, however, on where to begin.

The number of pulses is not known before hand thus, as SymPy is a symbolic solver, I must be able to dynamically define variables t1, t2, t3... tn.

My rough outline for an implementation is as follows:

import sympy as sp

def find_transmitter(arrival_times):
    equations = []

    for receiver in arrival_times:
        for pulse in arrival_times:
            equations.append(sp.Eq(c**2(pulse - ti...

    sp.solve(equations...

Where arrival_times is a list, where each column corresponds to a receiver, and each row corresponds to a pulse. I.e.:

[[a11, a21, a31...],[a12, a22, a32...],[a13, a23, a33...]]

It doesn't occur to me that there's a nice alternative to this approach. I'm not sure, however, e.g. how to dynamically define enough ti's and iterate through them. That is, e.g., I'll have an unknown number of equations looking like:

equations.append(sp.Eq(c**2(pulse - t1...
equations.append(sp.Eq(c**2(pulse - t2...
equations.append(sp.Eq(c**2(pulse - t3...
equations.append(sp.Eq(c**2(pulse - t4...
equations.append(sp.Eq(c**2(pulse - t5...
.
.
.

Is there a nice way to do this? Very likely, I'm overlooking some obvious solution.

10GeV
  • 453
  • 2
  • 14
  • If `arival_times` is large I'd expect very long `sympy.solve` times. `sympy` doesn't scale as well as numeric solvers (`numpy` and `scipy`). – hpaulj Nov 02 '21 at 21:50

2 Answers2

1

Perhaps surprisingly, can dynamically create any number of symbols with symbols using space as a delimiter!

You can also subscript with _N

arrival_times_count = 100  # from data length
chunksize = 10  # some factor of arrival_times_count

def chunks(lst, n):
    """Yield successive n-sized chunks from lst."""
    for i in range(0, len(lst), n):
        yield lst[i:i + n]

from sympy import *
# create a huge collection
collection = symbols(" ".join(f"t_{x}" for x in range(arrival_times_count)), real=True)
# fold it into rows
collection = list(chunks(collection, chunksize))

>>> for row in collection:
...     print(row)
...
(t_0, t_1, t_2, t_3, t_4, t_5, t_6, t_7, t_8, t_9)
(t_10, t_11, t_12, t_13, t_14, t_15, t_16, t_17, t_18, t_19)
(t_20, t_21, t_22, t_23, t_24, t_25, t_26, t_27, t_28, t_29)
(t_30, t_31, t_32, t_33, t_34, t_35, t_36, t_37, t_38, t_39)
(t_40, t_41, t_42, t_43, t_44, t_45, t_46, t_47, t_48, t_49)
(t_50, t_51, t_52, t_53, t_54, t_55, t_56, t_57, t_58, t_59)
(t_60, t_61, t_62, t_63, t_64, t_65, t_66, t_67, t_68, t_69)
(t_70, t_71, t_72, t_73, t_74, t_75, t_76, t_77, t_78, t_79)
(t_80, t_81, t_82, t_83, t_84, t_85, t_86, t_87, t_88, t_89)
(t_90, t_91, t_92, t_93, t_94, t_95, t_96, t_97, t_98, t_99)

chunks function from How do you split a list into evenly sized chunks?

However, there's a good chance this is an XY problem and you may be better off doing something like packing your data into a pandas.DataFrame, lambdifying your equation, and then calling .apply() to find the result of each row

ti7
  • 16,375
  • 6
  • 40
  • 68
1

A dirty way to dynamically define variables is to make a string pattern (with python code) and then process it with exec. So you can have access to the variables without passing though list, tuples, ... The best way would be to use substitution (vector -> components), hence you work with a vector & lambdify to perform componentwise operations.

from sympy import *
n = 4
t_s = ( ''.join((i, str(j))) for i, j in  zip('t'*n, range(1, n+1)))
exec( f"{', '.join(t_s)}  = symbols('t1:{n+1}')")

print(t2, type(t2))

Output

t2 <class 'sympy.core.symbol.Symbol'>
cards
  • 3,936
  • 1
  • 7
  • 25
  • honestly, dynamically creating sympy names/expressions is a great candidate for all sorts of eccentric `exec()` calls and really they go very well together and are something other languages don't have as easy! – ti7 Nov 03 '21 at 17:48
  • ...but you can fall also easily in your own trap if you are not careful! For example, `t0` will give raise error because the variables start at 1 and this is not very readable in the source code – cards Nov 03 '21 at 17:57
  • 1
    of course! still, you can inspect the state with _more evil_ like `set(globals().keys()) - previously_saved_keys` and `NameError` is easy to understand and handle! – ti7 Nov 03 '21 at 18:00