stiff.py

import numpy as np 
import matplotlib.pyplot as plt 

# 範囲
X = [0,1.5]

# 刻み幅
h1 = 0.5
h2 = 1.974/50
h3 = 1.875/50

x1 = np.arange(X[0],X[1]+h1,h1)
Nx1= len(x1)

x2 = np.arange(X[0],X[1]+h2,h2)
Nx2= len(x2)

x3 = np.arange(X[0],X[1]+h3,h3)
Nx3= len(x3)

# 導関数
f = lambda x,y : - 50*(y-np.cos(x))
# 初期値
y0 = 0.0

# Backward Euler
h = h1
Nx= Nx1
x = x1

y = np.zeros(Nx)
y[0] = y0

for n in range(Nx-1):
    y[n+1] = 1 / ( 1 + 50*h ) * (y[n] + 50 * h * np.cos(x[n+1]))
y_BE = y
x_BE = x

# Forward Euler1
h = h2
Nx= Nx2
x = x2

y = np.zeros(Nx)
y[0] = y0

h = h2
for n in range(Nx-1):
    y[n+1] = y[n] + h * f(x[n],y[n])
y_FE1 = y
x_FE1 = x

# Forward Euler2
h = h3
Nx= Nx3
x = x3

y = np.zeros(Nx)
y[0] = y0

h = h3
for n in range(Nx-1):
    y[n+1] = y[n] + h * f(x[n],y[n])
y_FE2 = y
x_FE2 = x

# plot
plt.figure()
plt.plot(x_FE1,y_FE1,'o-',linewidth=1.5,label='Forward Euler  (h=1.974/50)')
plt.plot(x_FE2,y_FE2,'o-',linewidth=1.5,label='Forward Euler  (h=1.875/50)')
plt.plot( x_BE, y_BE,'o-',linewidth=1.5,label='Backward Euler (h=0.5)')
plt.title('y\'=-50(y-cos(x))')
plt.legend(fontsize=12)
plt.grid()
plt.xlabel('x')
plt.ylabel('y')
plt.ylim([0,1.5])
plt.xlim([0,1.5])
plt.tight_layout()
plt.savefig('stiff.png')