drawio2py/experiments/Integrators/integrators_experiment.py

#!/usr/bin/python3
# This file was automatically generated from drawio2cbd with the command:
#   C:/Users/randy/Documents/git/DrawioConvert/__main__.py -F CBD -d experiments/Integrators -e G experiments/Integrators/integrators.drawio

from integrators_Integrators import *
from integrators_G import *
from integrators_Stiff import *
from CBD.simulator import Simulator

import matplotlib.pyplot as plt
import numpy as np

def func(t):
	return t / (t*t + 1)

def I(x):
	return np.log(x*x + 1) / 2

def stiff(t):
	return np.exp(-15 * t)

T = 1

times = np.arange(0, T, 0.001)
ys = stiff(times)

#plt.plot(times, ys, label="real")

# ints = {
# 	"forwards euler": [],
# 	"backwards euler": [],
# 	"trapezoid": [],
# 	"simpson 1/3": [],
# 	"analytic": 4.60522018
# }

from CBD.preprocessing.butcher import ButcherTableau as BT
from CBD.preprocessing.rungekutta import RKPreprocessor
cbd = Stiff("stiff")
cbd.addFixedRateClock("clock", 0.1)
tableau = BT.MidpointEuler()
RKP = RKPreprocessor(tableau, atol=2e-5, hmin=0.01, safety=.84)
transformed = RKP.preprocess(cbd)
sim = Simulator(transformed)
sim.run(T)

hist = transformed.getSignalHistory("y")
T = [t for t, _ in hist]
plt.plot(T, [x for _, x in hist], c="b")
plt.plot(T, [x for _, x in hist], 'o', fillstyle='none', lw=0.5, c="r")
plt.plot(T, [-0.1 for _ in range(len(hist))], '|', lw=0.5, c="g")
DT = [T[i + 1] - T[i] for i in range(len(T) - 1)]
print(min(DT), max(DT))

# for dt in [0.1, 0.01]:
# 	print("DT =", dt)
#
# 	cbd = G("G")
# 	# cbd = Stiff("stiff")
#
# 	# Run the Simulation
# 	sim = Simulator(cbd)
# 	sim.setDeltaT(dt)
# 	sim.run(T)
#
# # 	hist = cbd.getSignalHistory("y")
# # 	plt.plot([t for t, _ in hist], [x for _, x in hist], label="dt =" + str(dt))
# #
# 	ints["forwards euler"].append(cbd.getSignalHistory("fw")[-1].value)
# 	ints["backwards euler"].append(cbd.getSignalHistory("bk")[-1].value)
# 	ints["trapezoid"].append(cbd.getSignalHistory("trp")[-1].value)
# 	ints["simpson 1/3"].append(cbd.getSignalHistory("smp")[-1].value)
#
# 	hist = cbd.getSignalHistory("fw")
# 	plt.plot([t for t, _ in hist], [x for _, x in hist], '+', label="forwards, dt=" + str(dt))
# # 	# hist = cbd.getSignalHistory("bk")
# # 	# plt.plot([t for t, _ in hist], [x for _, x in hist], '.', label="backwards")
# # 	# hist = cbd.getSignalHistory("trp")
# # 	# plt.plot([t for t, _ in hist], [x for _, x in hist], 'x', label="trapezoid")
# 	hist = cbd.getSignalHistory("smp")
# 	plt.plot([t for t, _ in hist], [x for _, x in hist], ',', label="simpson, dt=" + str(dt))
# # 	# plt.plot(times, ys, label="real")
# # 	# plt.legend()
# # 	# plt.show()


# print("Results:")
# for k, v in ints.items():
# 	print(k + ":")
# 	print("\t", v)
#
# plt.legend()
plt.show()