List comprehension multiplication not expected

Question

from math import sin

def euler(f, x0, t0, h, N):
   t = t0
   x = x0
   while t <= N:
      t += h
      x += [h * x  for x in f(t, x)]
      print(x)

def f(t, x):
   vv = [-x[0]**3 - x[0] + sin(t)]
   return vv

This is my code. f is a function, x0 is the initial condition at time t0, t0 is the initial time, his the stepsize, and N is the number of steps. When I enter >>>euler(f, [0.], 0., 1., 10) I get [0.0, 0.8414709848078965, 0.9092974268256817, 0.1411200080598672, -0.7568024953079282, -0.9589242746631385, -0.27941549819892586, 0.6569865987187891, 0.9893582466233818, 0.4121184852417566, -0.5440211108893698] Which is incorrect.

I know my list comprehension statement is missing something, but I can't really pinpoint what I'm missing because I don't know how it is returning those values. The first 2 values are correct, then after that, it incorrectly calculates the remaining values.

When I am supposed to get

0.0
0.8414709848078965
0.313474190234726
0.11031613198378529
-0.758145003910546
-0.5231547727660838
-0.1362327890906342
0.6595149938422956
0.7024955870270317
0.06543687806493725

using this code:

from math import sin
def euler(f, x0, t0, h, N):
   t = t0
   x = x0
   while t <= N:
       t += h
       x += h * f(t,x)
       print(t,x)

def f(t, x):
    vv = -x**3 - x + sin(t)
    return vv

I think you mean to use the last value of `x` instead of the first, try replacing `x[0]` by `x[-1]` in `f()` — Julien, Nov 11 '16 at 09:39
I don't see what's the benefit of a listcomp here since you're iterating on 1 element. — Jean-François Fabre, Nov 11 '16 at 09:40
moreover, you're trying to use cumulative values ( relying on previous t and x) for which list comprehension may not be best suited. http://stackoverflow.com/questions/794774/python-list-comprehension-access-last-created-element explains it well — algrebe, Nov 11 '16 at 09:44
No it doesn't work. I get really high values, up to -1200. I'm more concerned with the list comprehension, because it does what it's supposed(oscillates) to but the values are off. — Anon Emouse, Nov 11 '16 at 09:46
@algrebe I'm not too advanced on this stuff, but from the page you linked me to, I'm guessing a for loop would be better? — Anon Emouse, Nov 11 '16 at 09:50
@AnonEmouse your current while loop is fine too. maybe you're looking for a way to make this more concise ? in which case you'd have better luck asking for a code review / alternate ways of doing it (i'm not sure whether stackoverflow is the right place though ). — algrebe, Nov 11 '16 at 09:56
No, I'm looking to fix the top code. It's not producing the right results. Because I want to be able to do it for multiple dimensions. I'm just not sure how to go about it.The bottom code works great but it's not exactly what I want. — Anon Emouse, Nov 11 '16 at 10:01
Oh ! Okay, I think you'll need to reword the question to say that you're trying to extend the bottom code to more than one dimension and provide an example of multidimensional input and output. — algrebe, Nov 11 '16 at 10:13
For truly multidimensional coupled ODE systems, use the `numpy` package. There is no need to re-invent a vector arithmetic library. -- Or if you really want to re-implement vector arithmetic based on the list type, to it properly and code a class `vector` for it, separating the vector internals from their application in the numerical code. — Lutz Lehmann, Nov 11 '16 at 10:19

score 2 · Accepted Answer · answered Nov 11 '16 at 10:15

Try this:

from math import sin

def euler(f, x0list, t0, h, N):
    t = t0
    xlist = x0list
    while t < N*h:
        klist = f(t,xlist)
        xlist = [x+h * k  for x,k in zip(xlist,klist)]
        t += h
        print t,xlist

def f(t, xlist):
    return [ -x**3 - x + sin(t) for x in xlist ]

euler(f, [0.,-0.1,0.1], 0., 1., 10)

This assumes that you want to evaluate multiple trajectories in parallel.

Note that the next step depends both on x and the slope k which again depends on x. Thus you also need a list evaluation in f.
With the zip pair building iterator the associated states and slopes are paired together for the computation of the next state.
Another way that avoids zip would have the function f return a list of pairs (x,k).

You could of course also avoid all this and iterate over the evaluation of euler where that only computes one trajectory.

For some reason, when I do `euler(f, [0.], 0., 0.2, 50)` it includes the `t` value of 10.19999 when in the code it clearly states `while t — Anon Emouse, Nov 11 '16 at 16:39
This is floating point.`0.2` is not representable as double, thus rounding errors occur and accumulate in `t+=h`. Replace `while ...` with `for _ in range(N):` to get an accurate number of loops. — Lutz Lehmann, Nov 11 '16 at 17:33

Lutz Lehmann · Answer 2 · 2016-11-11T11:00:30.660

If what you want is a function that returns the trajectory as a list, you can use the list generator pattern provided by yield:

def euler(f, x, t, h, N):
    yield x # report the initial point as the first sample
    for k in range(N):
        x += h*f(t,x)
        t += h
        yield x

def f(t, x):
    return -x**3 - x + sin(t)

x0, t0, tf = 0.0, 0.0, 10.0
N = 10
h = (tf-t0)/N

for x in euler(f, x0, t0, h, N):
     print x

for k,x in enumerate(euler(f, x0, t0, h, N)):
    print "%3d  %15.8f %20.16f" % (k, t0+k*h, x)

Making the method call more flexible in the form "tlist in - xlist out", use something like

def eulerfix(f, tlist, x):
    yield x
    for t, tnext in zip(tlist, tlist[1:]):
        x = x + (tnext-t)*f(t,x)
        yield x

def f(t, x):
    return -x**3 - x + sin(t)

x0, t0, tf = 0.0, 0.0, 10.0
N = 10
t = [ t0 + k*(tf-t0)/N for k in range(N+1) ]

for k,x in enumerate( eulerfix(f, t, x0) ):
    print "%3d  %15.8f %20.16f" % (k, t[k], x)

List comprehension multiplication not expected

2 Answers2