Search code examples
pythonnumerical-methodsode

Numerical stability of Euler's Method


I am trying to estimate when the following system of equations is numerically stable:

dx/dt = -x + y
dy/dt = -100*y

Solving for h in this inequality

|1+h*λ| <=1. where λ=-1,-100

(The eigenvalues of the system of equations). I get that for 1/50>= h>= 0. The Euler solution should be numerically stable. However, when I am plotting the numerical and analytical solution for timestep h=1/50. They don't seem to agree (while for h=1e-3 the results look good). Am I doing something wrong? This is the code I am using:

size=20
import matplotlib.pyplot as plt

# Define generl odel

w=1

def Euler(dt,y,x):
        
    dx = dt*(-x+y)
    dy     = -100*y*dt  
    y      = y+dy
    x  = x +dx
    return x,y


# Solve Diff Ep.

tstop =10.0
dt=1/50

y=np.zeros(int(tstop/dt))
x=np.zeros(int(tstop/dt))
t=np.zeros(int(tstop/dt))

t[0]     = 0
x[0]     = 0.1
y[0]     = 0.1



for i in range(0,int(tstop/dt)-1):
        
    dx     = dt*(-x[i]+y[i])
    dy     = -100*y[i]*dt  
    y[i+1]   = y[i]+dy
    x[i+1]   = x[i] +dx
    t[i+1]   = t[i]+dt
   

ax=plt.subplot(1,2,1)
plt.plot(t1,s[:,0],'g-', linewidth=3.0,label='X, Analytical solution')
plt.plot(t,x,"r--",label='X, Euler Method')
plt.plot(t1,s[:,1],'black',linewidth=3.0,label='Y, Analytical solution')
plt.plot(t,y,"b--",label='Y, Euler Method')
plt.text(0.2,0.1, r'$ {\Delta}t \ = \ $'+str(round(dt,4))+'s',horizontalalignment='center',
     verticalalignment='center', transform = ax.transAxes, fontsize=size)

ax.tick_params(axis='x',labelsize=17)
ax.tick_params(axis='y',labelsize=17)
plt.ylabel('X', fontsize=size)   
plt.xlabel('t [sec]', fontsize=size) 
plt.ylim([0,2])
plt.legend(frameon=False,loc=1,fontsize=18)
plt.xscale('log')


plt.ylim([0,0.12])

ax1.tick_params(axis='x',labelsize=17)
ax1.tick_params(axis='y',labelsize=17)
plt.rcParams["figure.figsize"] = [8,8]

enter image description here


Solution

  • With λ=-100 and h=1/50 you get for the propagation factor of the Euler method for the y component the value 1+h*λ=-1. Note that the solution, while oscillating, stays bounded. Which is the definition of A-stability. To get convergence to zero, you need h<0.02. To get non-oscillating behavior, you need h<=0.01. To see a somewhat curved behavior you need h<=0.005.