import numpy as np import matplotlib.pyplot as plt T=100*2*np.pi h = 2*np.pi/100 N=np.int(T/h) tr = np.linspace(0,T,N+1,endpoint=True) plt.figure(1) plt.subplot(1,2,1) plt.plot(tr,np.cos(tr)) plt.xlim([0,10*np.pi]) plt.subplot(1,2,2) plt.plot(tr,np.sin(tr)) plt.xlim([0,10*np.pi]) U= np.zeros((N+1,2)) U[0,0]=1.0 U[0,1]=0.0 for k in range(N): #Explicit Euler scheme U[k+1,0]=U[k,0] - h * U[k,1] U[k+1,1]=U[k,1] + h * U[k,0] plt.figure(2) plt.subplot(1,2,1) plt.title("EE real part") plt.plot(tr,np.cos(tr),"b-",tr,U[:,0],"r--") plt.xlim([0,10*np.pi]) plt.subplot(1,2,2) plt.title("EE imag part") plt.plot(tr,np.cos(tr),"b-",tr,U[:,1],"r--") plt.xlim([0,30*np.pi]) plt.ylim([-30,30]) U= np.zeros((N+1,2)) U[0,0]=1.0 U[0,1]=0.0 for k in range(N): #implicit Euler scheme U[k+1,0]=(U[k,0] - h * U[k,1])/(1+h*h) U[k+1,1]=U[k,1] + h * U[k+1,0] plt.figure(3) plt.subplot(1,2,1) plt.title("EI real part") plt.plot(tr,np.cos(tr),"b-",tr,U[:,0],"r--") plt.xlim([0,10*np.pi]) plt.subplot(1,2,2) plt.title("EI imag part") plt.plot(tr,np.cos(tr),"b-",tr,U[:,1],"r--") plt.xlim([0,30*np.pi]) plt.ylim([-30,30])