Forked from
DryVRgroup / DryVRtool
5 commits behind the upstream repository.
PDE.py 2.02 KiB
from scipy.integrate import odeint
import numpy as np
# import matplotlib.pyplot as plt
def PDE_dynamics(y,t,u1):
x1, x2, x3, x4, x5, x6 = y
x1 = float(x1)
x2 = float(x2)
x3 = float(x3)
x4 = float(x4)
x5 = float(x5)
x6 = float(x6)
x1_dot = 4.213*x3 - 252.57*x1 - 106.6*x2 - 51.94*u1 + 13.09*x4 - 1.8157*x5 + 0.4038*x6
x2_dot = 63.509*x3 - 106.6*x1 - 777.63*x2 - 11.124*u1 + 184.91*x4 - 26.094*x5 + 5.8051*x6
x3_dot = 21.751*x5 - 4.213*x1 - 63.509*x2 - 25.034*x3 - 251.84*x4 - 0.43288*u1 - 4.797*x6
x4_dot = 1.3463*u1 + 13.09*x1 + 184.91*x2 + 251.84*x3 - 634.37*x4 + 172.68*x5 - 39.239*x6
x5_dot = 172.68*x4 - 1.8157*x1 - 26.094*x2 - 21.751*x3 - 0.1867*u1 - 645.44*x5 + 337.53*x6
x6_dot = 39.239*x4 - 0.4038*x1 - 5.8051*x2 - 4.797*x3 - 0.041519*u1 - 337.53*x5 - 213.54*x6
dydt = [x1_dot, x2_dot, x3_dot, x4_dot, x5_dot, x6_dot]
return dydt
def TC_Simulate(Mode,initialCondition,time_bound):
time_step = 0.005;
time_bound = float(time_bound)
number_points = int(np.ceil(time_bound/time_step))
t = [i*time_step for i in range(0,number_points)]
if t[-1] != time_step:
t.append(time_bound)
newt = []
for step in t:
newt.append(float(format(step, '.3f')))
t = newt
u1 = 0.75
sol = odeint(PDE_dynamics, initialCondition, t, args=(u1,), hmax=time_step)
# Construct the final output
trace = []
for j in range(len(t)):
#print t[j], current_psi
tmp = []
tmp.append(t[j])
tmp.append(float(sol[j,0]))
tmp.append(float(sol[j,1]))
tmp.append(float(sol[j,2]))
tmp.append(float(sol[j,3]))
tmp.append(float(sol[j,4]))
tmp.append(float(sol[j,5]))
trace.append(tmp)
return trace
if __name__ == "__main__":
sol = TC_Simulate("Default", [-0.002, -0.0004, -0.001, -0.0019, -0.0008, -0.0001, 0.1, 1.0], 20.0)
#for s in sol:
# print(s)
time = [row[0] for row in sol]
a = [row[1] for row in sol]
plt.plot(time, a, "-r")
plt.show()
plt.show()