Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 57f5beba by Jigyasa Watwani
Options

deleting redundant files

parent e37a1c16
Inline Side-by-side
Showing with 0 additions and 414 deletions
linear_control_problem.py deleted 100644 → 0
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy.integrate import quad
from math import exp, sqrt, pi
def k_integral(l, l0, s, t):
return (np.exp((4 * l**2 + l0**2)/(8 * (s - t))) /(np.sqrt(np.pi) * 16 * (t - s)**(5/2)) *
( - ( (l0 - 2*l)**2 + 8 * (s - t)) * np.exp( (l0 + 2 * l)**2/(16 * (t - s)))
+ ( (l0 + 2*l)**2 + 8 * (s - t)) * np.exp( (l0 - 2 * l)**2/(16 * (t - s)))
)
)
def total_function(tt, l0, l):
ss = np.arange(0, tt-h/10, h)
index = np.arange(len(ss))
# print(l[index])
return np.trapz(k_integral(l[index], l0, ss[index], tt), ss[index])
# Explicit Euler Method
h = 1/(2**10)
t = np.arange(0, 5, h)
length = np.zeros(len(t))
l0 = 1
length[0] = l0
for i in range(0, len(t)-1):
length[i+1] = length[i] + h * total_function(t[i], l0, length) # length[i + 1] = length[i] + h * total_function(t[i], l0, length)
fig1, ax1 = plt.subplots(1, 1)
ax1.scatter(t, length)
ax1.set_ylabel("$L(t)$")
ax1.set_xlabel("$t$")
# # fig2, ax2 = plt.subplots(1, 1)
# #
# ax2.loglog(t, length)
# ax2.set_ylabel("$\log L(t)$")
# ax2.set_xlabel("$\log t$")
# fig3, ax3 = plt.subplots(1, 1)
# ax3.semilogy(t, length)
# ax3.set_ylabel("$\log L(t)$")
# ax3.set_xlabel("$t$")
plt.show()
print(length[-1])
\ No newline at end of file
linear_control_problem_an.py deleted 100644 → 0
import numpy as np
import matplotlib.pyplot as plt
import scipy
import math
import sympy
N = 40
l0 = 1
sigma0 = -1
def hermite_function(order, arg):
return ((2**order * math.factorial(order) * np.sqrt(np.pi))**(-0.5)
* scipy.special.eval_hermite(order, arg) * np.exp(-arg**2/2)
)
# times
T = 50
dt = 0.01
times = np.arange(0, T, dt)
# b
b = np.zeros((len(times), N), dtype='complex_')
b[0] = 0
# l
length = np.zeros(len(times))
length[0] = 1
# q1
q1_array = np.zeros(N, dtype='complex_')
q1_array[0] = ((1/2j) * ( (-1j) * np.sqrt(1/2) * (hermite_function(1, -2*l0) - hermite_function(1, 2*l0))
)
- l0 * ( (1/2) * (hermite_function(0, -2*l0) + hermite_function(0, 2*l0))
+ (-1j)**2 * np.sqrt(2)/2 * (hermite_function(2, -2*l0) + hermite_function(2, 2*l0))
)
)
q1_array[1] = ((1/2j) * ( np.sqrt(1/2) * (hermite_function(0, -2*l0) - hermite_function(0, 2*l0))
+ (-1j)**(2) * np.sqrt(2/2) * (hermite_function(2, -2*l0) - hermite_function(2, 2*l0))
)
- l0 * ( + (-1j) * (3/2) * (hermite_function(1, -2*l0) + hermite_function(1, 2*l0))
+ (-1j)**3 * np.sqrt(6)/2 * (hermite_function(3, -2*l0) + hermite_function(3, 2*l0))
)
)
for z in range(2, N):
q1_array[z] = ((1/2j) * ( (-1j)**(z-1) * np.sqrt(z/2) * (hermite_function(z-1, -2*l0) - hermite_function(z-1, 2*l0))
+ (-1j)**(z+1) * np.sqrt((z+1)/2) * (hermite_function(z+1, -2*l0) - hermite_function(z+1, 2*l0))
)
- l0 * ( (-1j)**(z-2) * np.sqrt((z-1)*z)/2 * (hermite_function(z-2, -2*l0) + hermite_function(z-2, 2*l0))
+ (-1j)**z * (2*z+1)/2 * (hermite_function(z, -2*l0) + hermite_function(z, 2*l0))
+ (-1j)**(z+2) * np.sqrt((z+1)*(z+2))/2 * (hermite_function(z+2, -2*l0) + hermite_function(z+2, 2*l0))
)
)
# q2
q2_array = np.zeros(N, dtype='complex_')
q2_array[0] =( (-1/4) * ((1/2) * (hermite_function(0, -3*l0/2) - hermite_function(0, -l0/2) - hermite_function(0, l0/2) + hermite_function(0, 3*l0/2))
+ (1/2) * np.sqrt(2) * (-1j)**(2) * (hermite_function(2, -3*l0/2) - hermite_function(2, -l0/2) - hermite_function(2, l0/2) + hermite_function(2, 3*l0/2))
)
- (l0/4j) * ((3/2) * np.sqrt(1/2) * (-1j) * (hermite_function(1, -3*l0/2) - hermite_function(1, -l0/2) - hermite_function(1, l0/2) + hermite_function(1, 3*l0/2))
+ (1/2) * np.sqrt(6) * (-1j)**3 * (hermite_function(3, -3*l0/2) - hermite_function(3, -l0/2) - hermite_function(3, l0/2) + hermite_function(3, 3*l0/2))
)
)
q2_array[1] = ( (-1/4) * ((3/2) * (-1j) * (hermite_function(1, -3*l0/2) - hermite_function(1, -l0/2) - hermite_function(1, l0/2) + hermite_function(1, 3*l0/2))
+ (1/2) * np.sqrt(6) * (-1j)**(3) * (hermite_function(3, -3*l0/2) - hermite_function(3, -l0/2) - hermite_function(3, l0/2) + hermite_function(3, 3*l0/2))
)
- (l0/4j) * ((3/2) * np.sqrt(1/2) * (hermite_function(0, -3*l0/2) - hermite_function(0, -l0/2) - hermite_function(0, l0/2) + hermite_function(0, 3*l0/2))
+ 3 * (-1j)**2 * (hermite_function(2, -3*l0/2) - hermite_function(2, -l0/2) - hermite_function(2, l0/2) + hermite_function(2, 3*l0/2))
+ (1/2) * np.sqrt(24) * (-1j)**4 * (hermite_function(4, -3*l0/2) - hermite_function(4, -l0/2) - hermite_function(4, l0/2) + hermite_function(4, 3*l0/2))
)
)
q2_array[2] = ( (-1/4) * ( (1/2) * np.sqrt(2) * (hermite_function(0, -3*l0/2) - hermite_function(0, -l0/2) - hermite_function(0, l0/2) + hermite_function(0, 3*l0/2))
+ (5/2) * (-1j)**2 * (hermite_function(2, -3*l0/2) - hermite_function(2, -l0/2) - hermite_function(2, l0/2) + hermite_function(2, 3*l0/2))
+ (1/2) * np.sqrt(12) * (-1j)**4 * (hermite_function(4, -3*l0/2) - hermite_function(4, -l0/2) - hermite_function(4, l0/2) + hermite_function(4, 3*l0/2))
)
- (l0/4j) * (3 * (-1j) * (hermite_function(1, -3*l0/2) - hermite_function(1, -l0/2) - hermite_function(1, l0/2) + hermite_function(1, 3*l0/2))
+ (9/2) * np.sqrt(3/2) * (-1j)**3 * (hermite_function(3, -3*l0/2) - hermite_function(3, -l0/2) - hermite_function(3, l0/2) + hermite_function(3, 3*l0/2))
+ (1/2) * np.sqrt(60) * (-1j)**5 * (hermite_function(5, -3*l0/2) - hermite_function(5, -l0/2) - hermite_function(5, l0/2) + hermite_function(5, 3*l0/2))
)
)
for zz in range(3, N):
q2_array[zz] = ((-1/4) * ((1/2) * np.sqrt(zz * (zz - 1)) * (-1j)**(zz - 2) * (hermite_function(zz-2, -3*l0/2) - hermite_function(zz-2, -l0/2) - hermite_function(zz-2, l0/2) + hermite_function(zz-2, 3*l0/2))
+ (2*zz + 1)/2 * (-1j)**zz * (hermite_function(zz, -3*l0/2) - hermite_function(zz, -l0/2) - hermite_function(zz, l0/2) + hermite_function(zz, 3*l0/2))
+ (1/2) * np.sqrt((zz + 1) * (zz + 2)) * (-1j)**(zz + 2) * (hermite_function(zz + 2, -3*l0/2) - hermite_function(zz + 2, -l0/2) - hermite_function(zz + 2, l0/2) + hermite_function(zz + 2, 3*l0/2))
)
- (l0 /4j) * ( (1/2) * np.sqrt(zz * (zz - 1) * (zz - 2)) * (-1j)**(zz - 3) * (hermite_function(zz - 3, -3*l0/2) - hermite_function(zz - 3, l0/2) - hermite_function(zz - 3, -l0/2) + hermite_function(zz - 3, 3*l0/2))
+ (3*zz / 2) * np.sqrt(zz/2) * (-1j)**(zz - 1) * (hermite_function(zz - 1, -3*l0/2) - hermite_function(zz - 1, l0/2) - hermite_function(zz - 1, -l0/2) + hermite_function(zz - 1, 3*l0/2))
+ (3*(zz + 1)/2) * np.sqrt((zz + 1) / 2) * (-1j)**(zz + 1) * (hermite_function(zz + 1, -3*l0/2) - hermite_function(zz + 1, l0/2) - hermite_function(zz + 1, -l0/2) + hermite_function(zz + 1, 3*l0/2))
+ (1/2) * np.sqrt((zz + 1) * (zz + 2) * (zz + 3)) * (-1j)**(zz + 3) * (hermite_function(zz + 3, -3*l0/2) - hermite_function(zz + 3, l0/2) - hermite_function(zz + 3, -l0/2) + hermite_function(zz + 3, 3*l0/2))
)
)
# q3
q3_array = np.zeros(N, dtype='complex_')
q3_array[0] = (1/4j) * (3*(1/2) * np.sqrt(1/2) * (-1j) * (hermite_function(1, 3*l0/2)
- hermite_function(1, l0/2)
+ hermite_function(1, -l0/2)
- hermite_function(1, -3*l0/2)
)
+ (1/2) * np.sqrt(3/2) * (-1j)**(3) * (hermite_function(3, 3*l0/2)
- hermite_function(3, l0/2)
+ hermite_function(3, -l0/2)
- hermite_function(3, -3*l0/2)
)
)
for yy in range(1, 3):
q3_array[yy] = (1/4j) * ( (3*yy/2) * np.sqrt(yy/2) * (-1j)**(yy-1) * (hermite_function(yy-1, 3*l0/2)
- hermite_function(yy-1, l0/2)
+ hermite_function(yy-1, -l0/2)
- hermite_function(yy-1, -3*l0/2)
)
+ (3*(yy+1)/2) * np.sqrt((yy+1)/2) * (-1j)**(yy+1) * (hermite_function(yy+1, 3*l0/2)
- hermite_function(yy+1, l0/2)
+ hermite_function(yy+1, -l0/2)
- hermite_function(yy+1, -3*l0/2)
)
+ (1/2) * np.sqrt((yy+1)*(yy+2)*(yy+3)/2) * (-1j)**(yy+3) * (hermite_function(yy+3, 3*l0/2)
- hermite_function(yy+3, l0/2)
+ hermite_function(yy+3, -l0/2)
- hermite_function(yy+3, -3*l0/2)
)
)
for y in range(3, N):
q3_array[y] = (1/4j) * ( (1/2) * np.sqrt((y*(y-1)*(y-2))/2) * (-1j)**(y-3) * (hermite_function(y-3, 3*l0/2)
- hermite_function(y-3, l0/2)
+ hermite_function(y-3, -l0/2)
- hermite_function(y-3, -3*l0/2)
)
+ (3*y/2) * np.sqrt(y/2) * (-1j)**(y-1) * (hermite_function(y-1, 3*l0/2)
- hermite_function(y-1, l0/2)
+ hermite_function(y-1, -l0/2)
- hermite_function(y-1, -3*l0/2)
)
+ (3*(y+1)/2) * np.sqrt((y+1)/2) * (-1j)**(y+1) * (hermite_function(y-1, 3*l0/2)
- hermite_function(y-1, l0/2)
+ hermite_function(y-1, -l0/2)
- hermite_function(y-1, -3*l0/2)
)
+ (1/2) * np.sqrt((y+1)*(y+2)*(y+3)/2) * (-1j)**(y+3) * (hermite_function(y+3, 3*l0/2)
- hermite_function(y+3, l0/2)
+ hermite_function(y+3, -l0/2)
- hermite_function(y+3, -3*l0/2)
)
)
# vector of hermite functions
vector_of_hermite_functions = np.zeros(N, dtype = 'complex_')
for j in range(0, N):
vector_of_hermite_functions[j] = hermite_function(j, 0) * (-1j)**(j)
# adjacency matrix
a = np.zeros((N, N), dtype='complex_')
for n in range(0, N):
for m in range(0, N):
if m == n:
a[m][n] = (n + 1/2)
elif m == n - 2:
a[m][n] = np.sqrt(n * (n - 1))/2
elif m == n + 2:
a[m][n] = np.sqrt((n + 1) * (n + 2))/2
# matrix
matrix = a + (1/np.pi) * np.dot(q1_array, vector_of_hermite_functions)
# T matrix
t = np.zeros((N+1, N+1), dtype='complex_')
for t1 in range(0, N):
t[N][t1] = (1/np.pi) * vector_of_hermite_functions[t1]
for t2 in range(0, N):
for t3 in range(0, N):
t[t2][t3] = matrix[t2][t3]
for t4 in range(0, N-1):
t[t4][N] = q3_array[t4]
# diagonalisation of t, p matrix
evals_t, evecs_t = np.linalg.eig(t)
diag_t = np.diag(evals_t)
d_inv = np.linalg.inv(diag_t)
p_inv = np.linalg.inv(evecs_t)
# exp of diag matrix
exp_diag_t = np.zeros((len(times), N+1, N+1), dtype='complex_')
for d1 in range(0, len(times)):
for nn in range(0, N):
exp_diag_t[d1][nn][nn] = np.exp(- times[d1] * evals_t[nn])
# q and q hat
q_array = np.append(q2_array, 1)
qhat = np.dot(p_inv, q_array)
# x
x = np.zeros((len(times), N+1) , dtype='complex_')
x[0] = np.append(b[0], length[0])
for tt in range(1, len(times)):
x[tt] = (np.linalg.multi_dot([evecs_t, exp_diag_t[tt], p_inv, x[0]])
+ 2 * sigma0 * (np.linalg.multi_dot([evecs_t, exp_diag_t[tt], d_inv, qhat])
- np.linalg.multi_dot([evecs_t, d_inv, qhat])
)
)
length[tt] = np.dot(np.append(np.zeros(N), 1), x[tt])
plt.plot(times, length)
plt.show()
\ No newline at end of file
linear_control_problem_hermite.py deleted 100644 → 0
import numpy as np
import matplotlib.pyplot as plt
import scipy
import math
N = 100
T = 5
Nt = 500
times = np.linspace(0, T, Nt)
dt = times[1] - times[0]
b = np.zeros((len(times), N))
b[0] = 0
q_array = np.zeros(N)
length = np.zeros(len(times))
length[0] = 1
l0 = length[0]
# matrix of hermite polynomials
matrix_of_hermite_polynomials = np.zeros(N, dtype = 'complex_')
for q in range(0, N):
normalization = (2**q * math.factorial(q) * np.sqrt(np.pi))**(-0.5)
matrix_of_hermite_polynomials[q] = normalization * scipy.special.eval_hermite(q, 0) * (-1j)**(q)
# adjacency matrix
a = np.zeros((N, N))
for n in range(0, N):
for m in range(0, N):
if m == n:
a[m][n] = (n + 1/2)
elif m == n - 2:
a[m][n] = np.sqrt(n * (n - 1))/2
elif m == n + 2:
a[m][n] = np.sqrt((n + 1) * (n + 2))/2
def q_array_function(length):
for p in range(0, N):
normalization = (2**p * math.factorial(p) * np.sqrt(np.pi))**(-0.5)
func = lambda k : 2 * (k**2 * np.sin(k * length) * np.sin(k * l0/2) * scipy.special.eval_hermite(p, k) * np.exp(-k**2/2))
q_array[p] = normalization * scipy.integrate.quad(func, -100, 100)[0]
return(q_array)
for i in range(0, len(times)-1):
b[i+1] = b[i] + dt * (q_array_function(length[i]) - np.dot(a, b[i]))
length[i+1] = length[i] + dt * (1/np.pi) * np.dot(matrix_of_hermite_polynomials, b[i])
# print(length[i])
plt.plot(times, length)
plt.show()
\ No newline at end of file
zeroth_control_problem_analytical.py deleted 100644 → 0
import numpy as np
import matplotlib.pyplot as plt
import scipy
import math
import sympy
N = 100
sigma0 = -1
l0 = 1
# t array
T = 100
dt = 0.01
times = np.arange(0, T, dt)
def hermite_function(order, arg):
return ((2**order * math.factorial(order) * np.sqrt(np.pi))**(-0.5)
* scipy.special.eval_hermite(order, arg) * np.exp(-arg**2/2)
)
# vector of hermite functions
vector_of_hermite_functions = np.zeros(N, dtype = 'complex_')
for q in range(0, N):
vector_of_hermite_functions[q] = hermite_function(q, 0) * (-1j)**(q)
# print(vector_of_hermite_functions)
# adjacency matrix
a = np.zeros((N, N), dtype='complex_')
for n in range(0, N):
for m in range(0, N):
if m == n:
a[m][n] = (n + 1/2)
elif m == n - 2:
a[m][n] = np.sqrt(n * (n - 1))/2
elif m == n + 2:
a[m][n] = np.sqrt((n + 1) * (n + 2))/2
# q1
q1_array = np.zeros(N, dtype='complex_')
q1_array[0] = (1/2j) * np.sqrt(1/2) * (-1j) * (hermite_function(1, -2 * l0) - hermite_function(1, 2 * l0))
for p1 in range(1, N):
q1_array[p1] = (1/2j) * ( np.sqrt(p1/2) * (-1j)**(p1-1) * (hermite_function(p1 - 1, -2 * l0) - hermite_function(p1 - 1, 2 * l0))
+ np.sqrt((p1 + 1)/2) * (-1j)**(p1+1) * (hermite_function(p1 + 1, -2 * l0) - hermite_function(p1 + 1, 2 * l0))
)
# q2
q2_array = np.zeros(N, dtype='complex_')
for p2 in range(0, 2):
q2_array[p2] = (-1/4) * ( ((2 * p2 + 1)/2) * (-1j)**p2 * (hermite_function(p2, -3*l0/2) - hermite_function(p2, -l0/2) - hermite_function(p2, l0/2) + hermite_function(p2, 3*l0/2))
+ np.sqrt((p2+1)*(p2+2))/2 * (-1j)**(p2+2) * (hermite_function(p2+2, -3*l0/2) - hermite_function(p2+2, -l0/2) - hermite_function(p2+2, l0/2) + hermite_function(p2+2, 3*l0/2))
)
for p3 in range(2, N):
q2_array[p3] = (-1/4) * ( ((2 * p3 + 1)/2) * (-1j)**p3 * (hermite_function(p3, -3*l0/2) - hermite_function(p3, -l0/2) - hermite_function(p3, l0/2) + hermite_function(p3, 3*l0/2))
+ np.sqrt((p3+1)*(p3+2))/2 * (-1j)**(p3+2) * (hermite_function(p3+2, -3*l0/2) - hermite_function(p3+2, -l0/2) - hermite_function(p3+2, l0/2) + hermite_function(p3+2, 3*l0/2))
+ np.sqrt(p3*(p3-1))/2 * (-1j)**(p3-2) * (hermite_function(p3-2, -3*l0/2) - hermite_function(p3-2, -l0/2) - hermite_function(p3-2, l0/2) + hermite_function(p3-2, 3*l0/2))
)
# matrix
matrix = a + (1/np.pi) * np.dot(q1_array, vector_of_hermite_functions)
# discard imaginary parts of matrix if they are < machine precision
index = np.where(np.imag(matrix) < 1e-10)
matrix[index] = np.real(matrix[index])
# print('M=', matrix)
# diagonal matrix
evalues_m, evectors_m = np.linalg.eig(matrix)
print(np.max(evalues_m), np.min(evalues_m))
diagonal_matrix = np.diag(evalues_m)
# print('D=', diagonal_matrix)
d_inv = np.linalg.inv(diagonal_matrix)
# print('$D^{-1}=$', d_inv)
q2_hat_array = np.dot(np.linalg.inv(evectors_m), q2_array)
# exp of diag matrix
exp_diag_matrix = np.zeros((len(times), N, N))
for tt in range(0, len(times)):
for nn in range(0, N):
# if (-times[tt] * evalues_m[nn]) < 2e01:
exp_diag_matrix[tt][nn][nn] = np.exp(- times[tt] * evalues_m[nn])
print(np.max(exp_diag_matrix), np.min(exp_diag_matrix))
# length
length = np.zeros(len(times))
for l in range(0, len(times)):
length[l] = l0 + (2 * sigma0/np.pi) * (np.linalg.multi_dot([vector_of_hermite_functions, evectors_m,
np.identity(N)*times[l], d_inv, q2_hat_array])
+ np.linalg.multi_dot([vector_of_hermite_functions, evectors_m, exp_diag_matrix[l],
d_inv, d_inv, q2_hat_array])
- np.linalg.multi_dot([vector_of_hermite_functions, evectors_m,
d_inv, d_inv, q2_hat_array])
)
# plt.plot(times, length)
# plt.show()
print(length[-1])
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment