length of the domain for delta function forcing at the centre of the domain -…

length of the domain for delta function forcing at the centre of the domain - hermite function way and integro-differential equation way for the linear problem

linear_control_problem.py

+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_hermite.py

+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