Initial Commit

Author: only-dev-ops

.gitignore

@@ -0,0 +1 @@
+.DS_Store

README.md

@@ -1,2 +1,9 @@
 # ODE-Methods
 A collection of various methods to find solution to Ordinary Differential Equations.
+
+1. Adams-Bashforth Method: Includes Two, Three, Four, Five Step method and also the Predictor-Corrector method
+2. Milne-Simspons Method
+3. Euler's Method: Includes the Analytical method, plain Estimation and Modified Euler's Method.
+4. Runge-Kutta Fourth-Order: Includes Estimator and Analytical method.
+5. Midpoint Method.
+6. Taylor's Estimation

ab_methods.py

@@ -0,0 +1,85 @@
+from math import *
+
+correctors = [
+	[2, 3, -1],
+	[12, 23, -16, 5],
+	[24, 55, -59, 37, -9],
+	[720, 1901, -2774, 2616, -1274, 251]
+]
+
+def abN_actual(f, f_actual, t, y, p, h, N):
+	corrector = correctors[N-2]
+	n = int((p-t)/h)
+	time = [t+h*i for i in range(n+1)]
+	values = [f_actual(time[i]) for i in range(N)]
+
+	for i in range(N, n+1):
+		fix = 0
+		for j in range(N):
+			fix += corrector[j+1] * f(time[i-j-1], values[i-j-1])
+		y = values[i-1] + h/corrector[0]*fix
+		values.append(y)
+
+	print(f"AB-{N} Method")
+	print("-"*40)
+	print(f't\t\t|\tEstimate')
+	print("-"*40)
+
+	for i in range(len(values)):
+		print(f'{time[i]:.3f}\t\t|\t{values[i]:.9f}')
+	print("-"*40)
+
+def abN_rk4(f, f_actual, t, y, p, h, N, correction = False):
+	corrector = correctors[N-2]
+	n = int((p-t)/h)
+	time = [t+h*i for i in range(n+1)]
+	y_values = [y]
+
+	for i in range(1,N):
+		y_values.append(rk4(f, time[i-1], y_values[i-1]))
+
+	for i in range(N, n+1):
+		fix = 0
+		for j in range(N):
+			fix += corrector[j+1] * f(time[i-j-1], y_values[i-j-1])
+		yn = y_values[i-1] + h/corrector[0]*fix
+
+		# This is a special case of AB-4 Predictor Corrector, hence the correction parameter for adjustment.
+		if N == 4 and correction:
+			fix = 9*f(time[i], yn)
+			coeffs = [19, -5, 1]
+			for j in range(len(coeffs)):
+				fix += coeffs[j] * f(time[i-j-1], y_values[i-j-1])
+			yn = y_values[i-1] + h/24*fix
+
+		y_values.append(yn)
+
+	print(f"AB-{N} {'Predictor-Corrector' if correction else ''} Method")
+	print("-"*70)
+	print(f't\t\t|\tEstimate\t\t|\tActual')
+	print("-"*70)
+
+	for i in range(len(y_values)):
+		print(f'{time[i]:.3f}\t\t|\t{y_values[i]:.9f}\t\t|\t{f_actual(time[i]):.9f}')
+	print("-"*70)
+
+def rk4(f, t, y):
+	k1 = h*f(t, y)
+	k2 = h*f(t + h/2, y + k1/2)
+	k3 = h*f(t + h/2, y + k2/2)
+	k4 = h*f(t + h, y + k3)
+	return y + (k1 + 2*k2 + 2*k3 + k4)/6
+
+if __name__ == '__main__':
+
+	f = lambda t, y: cos(2*t) + sin(3*t)
+
+	f_actual = lambda t: sin(2*t)/2 - cos(3*t)/3 + 4/3
+
+	t = float(input("t0: "))
+	y = float(input("y0: "))
+	p = float(input("Evaluation point: "))
+	h = float(input("Step size: "))
+	N = int(input("Which AB-N method to use? [2-5]: "))
+	#abN_actual(f, f_actual, t, y, p, h, N)
+	abN_rk4(f, f_actual, t, y, p, h, N, True)

eulers.py

@@ -0,0 +1,88 @@
+from math import *
+
+def analytical_eulers(f, f_actual, t, y, p, h = None, n = None):
+	if h == None:
+		h = (p-t)/n
+	if n == None:
+		n = int((p-t)/h)
+
+	print(f't\t\t|\tEstimate\t\t|\tExact\t\t\t|\tError')
+	print("-"*100)
+	print(f'{t:.3f}\t\t|\t{y:.9f}\t\t|\t{f_actual(t):.9f}\t\t|\t{abs(f_actual(t)-y):.9f}')
+
+	for i in range(n):
+		y = y + h*f(t, y)
+		t += h
+		actual = f_actual(t)
+		print(f'{t:.3f}\t\t|\t{y:.9f}\t\t|\t{actual:.9f}\t\t|\t{abs(actual-y):.9f}')
+
+	print("-"*100)
+
+def modified_eulers(f, f_actual, t, y, p, h = None, n = None):
+	if h == None:
+		h = (p-t)/n
+	if n == None:
+		n = int((p-t)/h)
+
+	print(f't\t\t|\tEstimate\t\t|\tExact\t\t\t|\tError')
+	print("-"*100)
+	print(f'{t:.3f}\t\t|\t{y:.9f}\t\t|\t{f_actual(t):.9f}\t\t|\t{abs(f_actual(t)-y):.9f}')
+
+	for i in range(n):
+		func_value = f(t, y)
+		t += h
+		next_value = f(t, y + h*func_value)
+		y = y + h/2*(func_value + next_value)
+		actual_value = f_actual(t)
+		print(f'{t:.3f}\t\t|\t{y:.9f}\t\t|\t{actual_value:.9f}\t\t|\t{abs(actual_value-y):.9f}')
+	print("-"*100)
+
+def eulers_estimate(f, t, y, p, h = None, n = None):
+	if h == None:
+		h = (p-t)/n
+	if n == None:
+		n = int((p-t)/h)
+
+	print(f't\t\t|\tEstimate')
+	print('-'*40)
+
+	print(f'{t:.3f}\t\t|\t{y:.9f}')
+
+	for i in range(n):
+		y = y + h*f(t, y)
+		t += h
+		print(f'{t:.3f}\t\t|\t{y:.9f}')
+	print('-'*40)
+
+def main():
+	option = int(input("Choose\n1) Euler's Approximation\n2) Euler's Analytical\n3) Modified Euler's\n"))
+
+	# Define these functions depending on the problem.
+	f = lambda t, y: -20*(y-t*t) + 2*t
+	f_analytical = lambda t: t*t + exp(-20*t)/3
+
+	t0 = float(input("t0 = "))
+	t = t0
+	y0 = float(input("y0 = "))
+	p = float(input("Evaluation point = "))
+	h = ""
+	n = ""
+	while h.strip() == "" and n.strip() == "":
+		h = input("Enter H [Leave blank if you want to enter N]: ")
+		n = input("Enter N [Leave blank if H already inputted]: ")
+	if h == "":
+		n = int(n)
+		h = None
+	else:
+		h = float(h)
+		n = None
+	if option == 1:
+		eulers_estimate(f, t0, y0, p, h, n)
+	elif option == 2:
+		analytical_eulers(f, f_analytical, t0, y0, p, h, n)
+	else:
+		modified_eulers(f, f_analytical, t0, y0, p, h, n)
+
+if __name__ == '__main__':
+	main()
+	

midpoint_method_ode.py

@@ -0,0 +1,43 @@
+from math import *
+
+def midpoint(f, f_actual, t, y, p, h = None, n = None):
+	if h == None:
+		h = (p-t)/n
+	if n == None:
+		n = int((p-t)/h)
+
+	print(f't\t\t|\tEstimate\t\t|\tExact\t\t\t|\tError')
+	print("-"*100)
+	print(f'{t:.3f}\t\t|\t{y:.9f}\t\t|\t{f_actual(t):.9f}\t\t|\t{abs(f_actual(t)-y):.9f}')
+
+	for i in range(n):
+		y = y + h*f(t + h/2, y + h/2*f(t, y))
+		t += h
+		actual = f_actual(t)
+		print(f'{t:.3f}\t\t|\t{y:.9f}\t\t|\t{actual:.9f}\t\t|\t{abs(actual-y):.9f}')
+	print("-"*100)
+
+def main():
+	# Define these functions depending on the problem.
+	f = lambda t, y: 1 + y/t
+	f_analytical = lambda t: t*(log(t) + 2)
+
+	t0 = float(input("t0 = "))
+	t = t0
+	y0 = float(input("y0 = "))
+	p = float(input("Evaluation point = "))
+	h = ""
+	n = ""
+	while h.strip() == "" and n.strip() == "":
+		h = input("Enter H [Leave blank if you want to enter N]: ")
+		n = input("Enter N [Leave blank if H already inputted]: ")
+	if h == "":
+		n = int(n)
+		h = None
+	else:
+		h = float(h)
+		n = None
+	midpoint(f, f_analytical, t0, y0, p, h, n)
+
+if __name__ == '__main__':
+	main()

milne_simpsons.py

@@ -0,0 +1,44 @@
+from math import *
+
+def milne_simpson(f, f_actual, t, y, p, h):
+	n = int((p-t)/h)
+	time = [t+i*h for i in range(n+1)]
+
+	y_values = [0]*(n+1)
+	y_values[0] = y
+
+	for i in range(1, 4):
+		y_values[i] = rk4(f, time[i-1], y_values[i-1])
+
+	for i in range(4, n+1):
+		milne = y_values[i-4] + 4/3*h*(2*f(time[i-1], y_values[i-1]) - f(time[i-2], y_values[i-2]) + 2*f(time[i-3], y_values[i-3]))
+		y_values[i] = y_values[i-2] + h/3*(f(time[i],milne) + 4*f(time[i-1], y_values[i-1]) + f(time[i-2], y_values[i-2]))
+
+	print(f"Milne-Simpson Method")
+	print("-"*70)
+	print(f't\t\t|\tEstimate\t\t|\tActual')
+	print("-"*70)
+
+	for i in range(len(y_values)):
+		print(f'{time[i]:.3f}\t\t|\t{y_values[i]:.9f}\t\t|\t{f_actual(time[i]):.9f}')
+	print("-"*70)
+
+def rk4(f, t, y):
+	k1 = h*f(t, y)
+	k2 = h*f(t + h/2, y + k1/2)
+	k3 = h*f(t + h/2, y + k2/2)
+	k4 = h*f(t + h, y + k3)
+	return y + (k1 + 2*k2 + 2*k3 + k4)/6
+
+if __name__ == '__main__':
+
+	f = lambda t, y: -5*y + 5*t*t + 2*t
+
+	f_actual = lambda t: t*t + exp(-5*t)/3
+
+	t = float(input("t0: "))
+	y = float(input("y0: "))
+	p = float(input("Evaluation point: "))
+	h = float(input("Step size: "))
+
+	milne_simpson(f, f_actual, t, y, p, h)

runge_kutta_4.py

@@ -0,0 +1,43 @@
+from math import *
+
+def rk4(f, f_actual, t, y, p, h, analytical = False):
+	n = int((p-t)/h)
+
+	separator = "-"*100 if analytical else "-"*40
+	if analytical:
+		print(f'\nt\t\t|\tEstimate\t\t|\tExact\t\t\t|\tError\n{separator}')
+		print(f'{t:.3f}\t\t|\t{y:.9f}\t\t|\t{f_actual(t):.9f}\t\t|\t{abs(f_actual(t)-y):.9f}')
+	else:
+		print(f'\nt\t\t|\tEstimate\n{separator}')
+		print(f'{t:.3f}\t\t|\t{y:.9f}')
+
+	for i in range(n):
+		k1 = h*f(t, y)
+		k2 = h*f(t + h/2, y + k1/2)
+		k3 = h*f(t + h/2, y + k2/2)
+		k4 = h*f(t + h, y + k3)
+		y = y + (k1 + 2*k2 + 2*k3 + k4)/6
+		t += h
+		actual = f_actual(t)
+		print(f'{t:.3f}\t\t|\t{y:.9f}', end = '')
+		print(f'\t\t|\t{actual:.9f}\t\t|\t{abs(actual-y):.9f}' if analytical else '')
+	print(separator)
+
+def main():
+
+	# Define these functions depending on the problem.
+	f = lambda t, y: -20*(y-t*t) + 2*t
+	f_analytical = lambda t: t*t + exp(-20*t)/3
+
+	choice = int(input("1) RK-4 Estimator\n2) RK-4 Analytical\n"))
+	if choice == 2:
+		print("Make sure you have defined the analytic function correctly.")
+
+	t0 = float(input("t0 = "))
+	y0 = float(input("y0 = "))
+	p = float(input("Evaluation point = "))
+	h = float(input("Step-Size = "))
+	rk4(f, f_analytical, t0, y0, p, h, choice == 2)
+
+if __name__ == '__main__':
+	main()

taylor_estimation.py

@@ -0,0 +1,27 @@
+from math import *
+
+g = 9.8
+kbym = 0.002/0.11
+print(kbym)
+
+yp = lambda t, y: -g - kbym*y*abs(y)
+ypp = lambda t, y, ypv: -kbym*(ypv*abs(y)+y*abs(ypv))
+
+series = lambda h, y, yp, ypp: y + h*yp + h*h*ypp/2
+
+a = 0
+b = 1
+t0 = 0
+y0 = 8
+h = 0.1
+n = int((b-t0)/h)
+
+for i in range(n):
+	d1 = yp(t0, y0)
+	d2 = ypp(t0, y0, d1)
+	print(f"y({t0:.3f}) = {y0}\ny'({t0:.3f}) = {d1}\ny''({t0:.3f}) = {d2}")
+	result = series(h, y0, d1, d2)
+	t0 += h
+	print(f'y({t0:.3f}) = {result:.9f}\n')
+	y0 = result
+