Runge-Kutta method

integrator.py

@@ -1,31 +1,82 @@
 from volume import eval_volume
 import matplotlib.pyplot as plt
+from enum import Enum
 
 
-def euler_integrator(system, values, t, h, use_modified_method=False):
+class Method(Enum):
+ EULER = 0
+ MODIFIED_EULER = 1
+ EULER_CROMER = 2
+ RUNGE_KUTTA_5 = 3
+
+
+def euler_step(system, values, h):
+ size = len(system)
+ temp_values = values.copy()
+ for i in range(size):
+ values[i] += system[i](temp_values) * h
+
+
+def modified_euler_step(system, values, h):
+ size = len(system)
+ temp_values = []
+ for i in range(size):
+ temp_values.append(values[i] + h / 2 * system[i](values))
+ for i in range(size):
+ values[i] += system[i](temp_values) * h
+
+
+def euler_cromer_step(system, values, h):
+ size = len(system)
+ for i in range(size):
+ values[i] += system[i](values) * h
+
+
+runge_kutta_a = [
+ [1/5],
+ [3/40, 9/40],
+ [44/45, -56/15, 32/9],
+ [19372/6561, -25360/2187, 64448/6561, -212/729],
+ [9017/3168, -355/33, 46732/5247, 49/176, -5163/18656],
+ [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]
+]
+
+
+def runge_kutta_5_step(system, values, h):
+ size = len(system)
+ k = []
+ temp_values = [values.copy()]
+ for rk_order in range(5):
+ k.append([])
+ temp_values.append(values.copy())
+ for i in range(size):
+ k[rk_order].append(system[i](temp_values[rk_order]))
+ for rk_line in range(rk_order + 1):
+ for i in range(size):
+ temp_values[rk_order + 1][i] += h * runge_kutta_a[rk_order][rk_line] * k[rk_order][i]
+ values[:] = temp_values[-1][:]
+
+
+def integrator(system, values, t, h, method):
  size = len(system)
  time_history = []
  values_history = []
  for i in range(size):
  values_history.append([])
  for step in range(0, int(t / h)):
  time_history.append(step * h)
- if use_modified_method:
- temp_values = []
- for i in range(size):
- temp_values.append(values[i] + h / 2 * system[i](values))
- for i in range(size):
- values[i] += system[i](temp_values) * h
- values_history[i].append(values[i])
- else:
- for i in range(size):
- values[i] += system[i](values) * h
- values_history[i].append(values[i])
+ match method:
+ case Method.EULER: euler_step(system, values, h)
+ case Method.MODIFIED_EULER: modified_euler_step(system, values, h)
+ case Method.EULER_CROMER: euler_cromer_step(system, values, h)
+ case Method.RUNGE_KUTTA_5: runge_kutta_5_step(system, values, h)
+ for i in range(size):
+ values_history[i].append(values[i])
  return time_history, values_history
 
 
-def run(system, initial_values, t, h, window, use_modified_method=False):
- time_history, values_history = euler_integrator(system, initial_values, t, h, use_modified_method)
+def run(system, initial_values, t, h, window, method):
+ time_history, values_history = integrator(system, initial_values, t, h, method)
  steps_history, volume_history = eval_volume(values_history, window)
  size = len(initial_values)
 

main.py

@@ -1,4 +1,4 @@
-from integrator import run
+from integrator import run, Method
 from gui import *
 
 
@@ -7,16 +7,18 @@
 window_text, window_field = create_field("Window")
 
 tab_group = SGui.TabGroup([
- [SGui.Tab("Default", create_layout(Layout.Default), key=Layout.Default)],
+ [SGui.Tab("Van der Pol", create_layout(Layout.Default), key=Layout.Default)],
  [SGui.Tab("Rossler", create_layout(Layout.Rossler), key=Layout.Rossler)]
  ], expand_x=True)
 
-modified_method_checkbox = SGui.Checkbox("Use Modified Method")
+methods = ["Euler", "Modified Euler", "Euler-Cromer", "Runge-Kutta 5"]
+method_key = "Method"
+method_field = SGui.Combo(methods, key=method_key)
 
 layout = [
  [t_text, t_field, h_text, h_field],
  [tab_group],
- [modified_method_checkbox],
+ [method_field],
  [window_text, window_field, SGui.OK("Run", expand_x=True)]
 ]
 app = SGui.Window("Numerical Integrator", layout)
@@ -30,6 +32,7 @@
  t, h, window = get_float(values, 'T'), get_float(values, 'h'), get_int(values, 'Window')
  x, y = get_float(values, 'x_' + layout_name), get_float(values, 'y_' + layout_name)
  app.hide()
+ method = Method(methods.index(values[method_key]))
  if tab_group.Get() == Layout.Default:
  mu = get_float(values, 'μ')
  run(
@@ -38,8 +41,7 @@
  lambda var: mu * (1 - pow(var[0], 2)) * var[1] - var[0] # mu * (1 - pow(x, 2)) * y - x
  ],
  [x, y],
- t, h, window,
- modified_method_checkbox.Get()
+ t, h, window, method
  )
  else:
  a, b, c = get_float(values, 'a'), get_float(values, 'b'), get_float(values, 'c')
@@ -50,8 +52,7 @@
  lambda var: b + var[2] * (var[0] - c) # b + z * (x - c)
  ],
  [x, y, get_float(values, 'z')],
- t, h, window,
- modified_method_checkbox.Get()
+ t, h, window, method
  )
  app.un_hide()