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]
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
.