Euler's method.

Question

So for my assignment I have to code a program to solve first order ODE's using Euler's Method. My program works, it returns the right values. (I checked using an online calculator). However, solution to the assignment returns something very different.

The initial condition is x(0) = 0, from t = 0 to t = 10. With 10 steps. (I presume it to be the green line)

Solution to answer

the function is this

Using my program I get the results: ([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], (the x values)

And the y-values respectively

[[0.0]

[-2.0]

[-12.909297426825681]

[-43.245542352772254]

[-110.57365446753619]

[-239.994320283207]

[-462.9378998529216]

[-812.0361590626763]

[-1333.0340119617658]

[-2071.8747943733247]

[-3080.8748098338056]

So as you can see, in the solutions, the y values don't even go as far as -1.0. So I'm wondering how does this happen?

EDIT CODE: f = function y0 = initial condition at time t0 t0 = initial time h = step size N = number of steps

def integrateEuler(f,y0,t0,h,N):
t = t0
y = y0
z = []
v = []
yf = N*h #final xval
while t <= N:
    xval = t
    yval = [y]
    t += h
    y += h * f(t,y)
    z.append(xval)
    v.append(yval)
return z, v
#FUNCTION
def f(x,t):
    vv = -x**3 - x + sin(t)
    return vv

I enter in the shell:

 >>>Euler(f, 0., 0., 1., 10)

Answer 1

This may be a familiar problem with forward Euler: it is not stable for arbitrarily large step sizes. It can be seen in the simpler problem $y'=-y,y(0)=1$. In this case, if $h>2$, then after one step forward in time, the numerical solution has not only changed sign (which it shouldn't be able to do), it has also increased in magnitude (which it also shouldn't be able to do). As a result, the numerical solution oscillates without bound as you move forward in time, while the true solution is $y=e^{-x}$, which of course decays as you move forward in time.

In your equation essentially the same thing is going on: the $-x^3-x$ part is trying to pull the true ODE back toward $0$ (you can see this by considering $\frac{d}{dt} \left ( \frac{x^2}{2} \right ) = -x^4-x^2+x \sin(t)$, which is negative for, say, $|x|>1$). But if the step size is too large then the numerical solution goes well past $0$ (when started at a positive value, say), which allows for unbounded oscillation to begin just like in the previous example. Because your equation also has a forcing, this problem can happen even when you begin at $x=0$.

It could also just be a programming error; I didn't check over the code.

Answer 2

Compare the lines

        y += h * f(t,y)

and

def f(x,t):

and after correcting this inconsistency (and some others relating to the choice of t value used), your values should track the ones in the plot.