CODEWITHSUNDEEP
pythonmatplotlibdifferential-equations
In the blind
In the blind

Reputation: 113

Euler's method for Python

I want to approximate the solutions of dy/dx = -x +1, with eulers method on the interval from 0 to 2. I'm using this code

def f(x):
    return -x+1  # insert any function here


x0 = 1  # Initial slope #
dt = 0.1  # time step
T = 2  # ...from...to T
t = np.linspace(0, T, int(T/dt) + 1)  # divide the interval from 0 to 2 by dt
x = np.zeros(len(t))
x[0] = x0 # value at 1 is the initial slope
for i in range(1, len(t)):  # apply euler method
    x[i] = x[i-1] + f(x[i-1])*dt

plt.figure()  # plot the result
plt.plot(t,x,color='blue')
plt.xlabel('t')
plt.ylabel('y(t)')


plt.show()

Can I use this code to approximate the solutions of any function on any interval? It's hard to see whether this actually works, because I don't know how to plot the actual solution ( -1/2x^2 + x ) along side the approximation.

Upvotes: 0

Views: 2176

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

It would probably help if you consistently used the same variable names for the same role. Per your output, the solution is y(t). Thus your differential equation should be dy(t)/dt = f(t,y(t)). This would then give an implementation for the slope function and its exact solution

def f(t,y): return 1-t

def exact_y(t,t0,y0): return y0+0.5*(1-t0)**2-0.5*(1-t)**2

Then implement the Euler loop also as a separate function, keeping out problem specific details as much as possible

def Eulerint(f,t0,y0,te,dt):
    t = np.arange(t0,te+dt,dt)
    y = np.zeros(len(t))
    y[0] = y0 
    for i in range(1, len(t)):  # apply euler method
        y[i] = y[i-1] + f(t[i-1],y[i-1])*dt
    return t,y

Then plot the solutions as

y0,T,dt = 1,2,0.1
t,y = Eulerint(f,0,y0,T,dt)
plt.plot(t,y,color='blue')
plt.plot(t,exact_y(t,0,y0),color='orange')

enter image description here

Upvotes: 2

ernesi
ernesi

Reputation: 378

You can just plot the actual solution by using:

def F(x):
   return -0.5*x+x

# some code lines

plt.plot(t,x,color='blue')
plt.plot(t,F(t),color='orange')

But please note that the actual solution (-1/2x+x = 1/2x) does not correspond to your slope f(x) and will show a different solution.

The *real slope f(x) of the actual solution (-1/2x+x = 1/2x) is just simply f(x)=1/2

Upvotes: 1

Related Questions