impulse forcing length vs time full non-linear problem

full_nonlinear_problem.py

+import numpy as np
+import matplotlib.pyplot as plt
+import scipy
+import math
+import sympy
+
+N = 80
+sigma0 = -1
+l0 = 1
+
+# t array
+T = 150
+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_polynomials = np.zeros(N, dtype = 'complex_')
+for q in range(0, N):
+    vector_of_hermite_polynomials[q] = hermite_function(q, 0) * (-1j)**(q) 
+
+# 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
+
+# a1
+a1 = np.zeros((N, N), dtype='complex_')
+for nn in range(0, N):
+    for mm in range(0, N):
+        if mm == nn-1:
+            a1[nn][mm] = np.sqrt(nn/2)
+        if mm == nn + 1:
+            a1[nn][mm] = np.sqrt((nn + 1)/2)
+
+# q1
+def q1(length):
+    q1_array = np.zeros(N, dtype='complex_')
+    q1_array[0] = np.sqrt(1/2) * (-1j) * (hermite_function(1, 0) + hermite_function(1, -2 * length))
+    for p1 in range(1, N):
+        q1_array[p1] = ( np.sqrt(p1/2) * (-1j)**(p1-1) * (hermite_function(p1 - 1, 0) + hermite_function(p1 - 1, -2 * length))
+                        + np.sqrt((p1 + 1)/2) * (-1j)**(p1+1) * (hermite_function(p1 + 1, 0) + hermite_function(p1 + 1, -2 * length))    
+                                )
+    return q1_array
+
+# q2
+def q2(length):
+    q2_array = np.zeros(N, dtype='complex_')
+    for p2 in range(0, 1):
+        q2_array[p2] = (1/2j) * (     ((2 * p2 + 1)/2)         * (-1j)**p2 *     (hermite_function(p2, length + l0/2) - hermite_function(p2, length-l0/2))
+                                    + np.sqrt((p2+1)*(p2+2))/2 * (-1j)**(p2+2) * (hermite_function(p2+2, length + l0/2) - hermite_function(p2+2, length - l0/2))
+                                )
+    for p3 in range(2, N):
+        q2_array[p3] = (1/2j) * (   ((2 * p3 + 1)/2)           * (-1j)**p3     * (hermite_function(p3, length+l0/2) - hermite_function(p3, length-l0/2))
+                                    + np.sqrt((p3+1)*(p3+2))/2 * (-1j)**(p3+2) * (hermite_function(p3+2, length+l0/2) - hermite_function(p3+2, length-l0/2))
+                                    + np.sqrt(p3*(p3-1))/2     * (-1j)**(p3-2) * (hermite_function(p3-2, length+l0/2) - hermite_function(p3-2, length-l0/2))
+                                )
+    return q2_array
+# b 
+b = np.zeros((len(times), N), dtype='complex_')
+b[0] = 0
+
+length = np.zeros(len(times))
+length[0] = l0
+
+def func1(b, length):
+    return ((1/np.pi) * np.dot(vector_of_hermite_polynomials, b) * (1j * q1(length) + np.dot(a1, b))
+            - np.dot(a, b) - 2j * sigma0 * q2(length)
+           )
+def func2(b):
+    return (1/np.pi) * np.dot(vector_of_hermite_polynomials, b)
+
+
+# rk4
+for l in range(0, len(times)-1):
+    b[l+1] = b[l] + dt * func1(b[l], length[l])
+    length[l+1] = length[l] + dt * func2(b[l])
+
+print(length[-1])
+plt.plot(times, length)
+plt.show()
\ No newline at end of file