length vs time exact solutions for impulse forcing at zeroth and linear order.…

length vs time exact solutions for impulse forcing at zeroth and linear order. Hasn't converged yet.

zeroth_control_problem_analytical.py

+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