|
# Runge-Kutta Method
import numpy as np import matplotlib.pyplot as plt
def RK2(f, h=1./16, y0=0, a=0, b=1): t = np.arange(a,b+h,h); y = np.zeros_like(t) y[0] = y0; N = t.size for n in range(0,N-1): v1 = f( t[n],y[n] ); v2 = f( t[n]+h, y[n]+h*v1 ) y[n+1] = y[n] + h/2*( v1 + v2 ) EY = Y(t) err = h**(1/2)*np.linalg.norm(EY-y) plt.plot(t,EY, t,y,'r--'); plt.pause(.5) print(' RK2-err = ', '%12.6e' %err, ' with h = ', h) return err def RK4(f, h=1/16, y0=0, a=0, b=1): t = np.arange(a,b+h,h); y = np.zeros_like(t) y[0] = y0; N = t.size for n in range(0,N-1): v1 = f( t[n],y[n] ) v2 = f( t[n]+h/2, y[n]+h/2*v1 ) v3 = f( t[n]+h/2, y[n]+h/2*v2 ) v4 = f( t[n]+h, y[n]+h*v3 ) y[n+1] = y[n] + h/6*( v1 + 2*v2 + 2*v3 + v4 ) EY = Y(t) err = h**(1/2)*np.linalg.norm(EY-y) plt.plot(t,EY, t,y,'r--'); plt.pause(.5) print(' RK4-err = ', '%12.6e' %err, ' with h = ', h) return err
def f(t,y): return -y + 2*np.cos(t) # return y/4*( 1- y/20 )
def Y(t): return np.sin(t) + np.cos(t) # return 20/( 1+19*np.exp(-t/4) )
# Solve y0 = 1; a= 0; b = np.pi; #y0 = 1; a= 0; b = 30; plt.clf() e1 = []; e2 = [] for k in range(1,10): h = 2**(-k) e1.append( RK2(f,h,y0,a,b) )
e1 = np.array(e1) print('RK2 °¨Ãà ºñÀ² = ', e1[:-1]/e1[1:])
for k in range(1,10): h = 2**(-k) e2.append( RK4(f,h,y0,a,b) ) e2 = np.array(e2) print('RK4 °¨Ãà ºñÀ² = ', e2[:-1]/e2[1:])
|
|
|
|
|
|