ODE
7.3) Write a fourth-order Runge–Kutta adaptive stepper for the preceding problem, and check how the average step size that it finds depends on the desired local error.
In [0]:
import math
from matplotlib import pyplot as plt
import numpy as np
In [0]:
history_x = []
history_h = []
history_y_rk4 = []
history_error_rk4 = []
#initial condition
current_y = 1
dy = 0
x = 0
h = 0.001
In [0]:
def y(t):
y = math.cos(t)
return y
def dy(t):
dy = -1*math.sin(t)
return dy
def ddy(t):
ddy = -1*math.cos(t)
return ddy
In [0]:
def next_y_rk4(t,current,h):
k1 = h*dy(t)
k2 = h*dy(t+0.5*h)
k3 = h*dy(t+0.5*h)
k4 = h*dy(t)
y_next_temp = current + k1/6 + k2/3 + k3/3 + k4/6
return y_next_temp
In [0]:
#step_Count = 500
max_Value = 10*math.pi
error_threshold = 0.0001
h_upper_test = 0
while(x<max_Value):
# Check Error
#print("---------loop start---------")
y_temp_1 = next_y_rk4(x,current_y,h)
y_temp_2 = next_y_rk4(x,current_y,0.5*h)
y_temp_2 = next_y_rk4(x+0.5*h,current_y, 0.5*h)
error_temp = abs(y_temp_2-y_temp_1)
#print("h is " + str(h))
#print("error is " + str(error_temp))
#print("Flag is " + str(h_upper_test))
if(error_temp < error_threshold):
# check bigger step h
if(h_upper_test == 1):
# this is right h just do normal h step
current_y = y_temp_1
x = x + h
h_upper_test = 0;
#print("go000000000000000000000000000000")
history_h.append(h)
history_x.append(x)
history_y_rk4.append(current_y)
else:
## I havent tried bigger h yet, lets try
h = h * 2.0
h_upper_test = 0;
#("h x 2 ")
else:
# current h is not good. change h to h*0.5
h = h * 0.5
h_upper_test = 1;
#print("h / 2 ")
# End While
In [6]:
tempCount = range(1,len(history_h)+1)
print("error threshold is " + str(error_threshold))
plt.plot(history_x,history_h,'b',c='red',label="h")
#plt.plot(history_h,history_error_rk4,'b',c='blue',label="rk4 error")
plt.xlabel('x')
plt.ylabel('h')
plt.legend()
plt.show()
plt.plot(history_x,history_y_rk4,'b',c='red',label="y")
#plt.plot(history_h,history_error_rk4,'b',c='blue',label="rk4 error")
plt.xlabel('x')
plt.ylabel('y_rk4')
plt.legend()
plt.show()
