I have a coupled pair of anharmonic oscillators and need to calculate the stability matrix and find the lyaupanov exponent for how the nearby trajectories diverge. In particular, I have the Hamiltonian
H = (p1^2+p2^2 + q1^4+q2^4 + 12*q1^2*q2^2 )/2
and I need to compute the matrix M given by
dM/dt = J*Hess*M
where J := Matrix(4,4,[0,0,1,0, 0,0,0,1, -1,0,0,0, 0,-1,0,0]), and the Hessian takes the form:
Matrix(4,4, [6*q1_12(t)^2 + 12*Q2_12(t)^2, 2*12*Q1_12(t)*Q2_12(t),0,0, 2*12*Q1_12(t)*Q2_12(t), 6*Q2_12(t)^2 + 12*q1_12(t)^2,0,0, 0,0,1,0, 0,0,0,1])
The solver finds the trajectories of p1,p2,q1,q2 fine, but I don't seem to find a way to incorporate their solutions as inputs to reevaluate the Hessian at each time step. I read through https://www.maplesoft.com/support/help/Maple/view.aspx?path=DEtools%2fmatrixDE and the dsolve since I'm trying to do this numerically.
I thought I could get around this just by resolving the trajectories, but it's spitting out an error that arrays must be initialized with lists.
The actual code is here: AMO_HW4.mw
Sorry, I still need to clean it up a bit. Any help would be appreciated. I need to calculate M(t), then calculate the matrix norm and find the exponent.