import numpy as np import matplotlib.pyplot as plt Y0 = 0.5 T0 = 0. T1 = 4. # Question 1 def f(y, t): return y-t*t+1 # Question 2 def RK4(f, Y0, t0, t1, h): T = [t0] Y = [Y0] y = Y0 t = t0 while t < t1: k1 = h*f(y, t) k2 = h*f(y + k1 / 2., t + h / 2.) k3 = h*f(y + k2 / 2., t + h / 2.) k4 = h*f(y + k3, t + h) y += 1/6.*(k1 + 2 * k2 + 2 * k3 + k4) t += h Y.append(y) T.append(t) return T, Y # Question 3 def f2(t): return t*t+2*t+1-0.5*np.exp(t) # Question 4 H = [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3] for h in H: T, Y = RK4(f, Y0, T0, T1, h) plt.plot(T, Y) T = np.arange(T0, T1, 0.1) plt.plot(T, f2(T), '*') plt.savefig('RK4.png')