JacCalculating Jacobian Matrix for System of Differential Equations

[bangthatdrum]

Greetings all,

I hope someone out there in the vast hinterland of the internet can help.

I'm trying to calculate the lyapunov exponent for a system of differential equations. Now I can do this just fine for a system involving only first order derivatives such the Lorenz system, however, and this is the kicker, I'm trying to calculate it for a system involving 2nd and 1st order derivatives.

So my question, what would the jacobian be for a systems like this?

\begin{align*}
\ddot{y}_1+\alpha\dot{y}_1+\beta y_1&=\beta (w_{ee}S(y_2) -w_{ie}S(y_3) -w_{ie}S(y_4) )\\
\ddot{y}_2+\alpha\dot{y}_2+\beta y_2&=\beta (w_{ee}S(y_1) -w_{ie}S(y_3))\\
\ddot{y}_3+\alpha\dot{y}_3+\beta y_3&=\beta (w_{ei}S(y_1) +w_{ei}S(y_2) -w_{ii}S(y_4))\\
\ddot{y}_4+\alpha\dot{y}_4+\beta y_4&=\beta (w_{ei}S(y_1) -w_{ii}S(y_3))
\end{align*}

I initially thought you just rearranged it for y1ddot and then (in Matlab code) did something like this, but now I'm not so sure this is the right arrangement. I've looked everywhere I can think of but can't find any info.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a11=diff(y1dot,y1);
a12=diff(y1dot,y1dot);
a13=diff(y1dot,y2);
a14=diff(y1dot,y2dot);
a15=diff(y1dot,y3);
a16=diff(y1dot,y3dot);
a17=diff(y1dot,y4);
a18=diff(y1dot,y4dot);

a21=diff(y1ddot,y1);
a22=diff(y1ddot,y1dot);
a23=diff(y1ddot,y2);
a24=diff(y1ddot,y2dot);
a25=diff(y1ddot,y3);
a26=diff(y1ddot,y3dot);
a27=diff(y1ddot,y4);
a28=diff(y1ddot,y4dot);

%--------------------------------------------------------------------------

a31=diff(y2dot,y1);
a32=diff(y2dot,y1dot);
a33=diff(y2dot,y2);
a34=diff(y2dot,y2dot);
a35=diff(y2dot,y3);
a36=diff(y2dot,y3dot);
a37=diff(y2dot,y4);
a38=diff(y2dot,y4dot);

a41=diff(y2ddot,y1);
a42=diff(y2ddot,y1dot);
a43=diff(y2ddot,y2);
a44=diff(y2ddot,y2dot);
a45=diff(y2ddot,y3);
a46=diff(y2ddot,y3dot);
a47=diff(y2ddot,y4);
a48=diff(y2ddot,y4dot);

%--------------------------------------------------------------------------

a51=diff(y3dot,y1);
a52=diff(y3dot,y1dot);
a53=diff(y3dot,y2);
a54=diff(y3dot,y2dot);
a55=diff(y3dot,y3);
a56=diff(y3dot,y3dot);
a57=diff(y3dot,y4);
a58=diff(y3dot,y4dot);

a61=diff(y3ddot,y1);
a62=diff(y3ddot,y1dot);
a63=diff(y3ddot,y2);
a64=diff(y3ddot,y2dot);
a65=diff(y3ddot,y3);
a66=diff(y3ddot,y3dot);
a67=diff(y3ddot,y4);
a68=diff(y3ddot,y4dot);

%--------------------------------------------------------------------------

a71=diff(y4dot,y1);
a72=diff(y4dot,y1dot);
a73=diff(y4dot,y2);
a74=diff(y4dot,y2dot);
a75=diff(y4dot,y3);
a76=diff(y4dot,y3dot);
a77=diff(y4dot,y4);
a78=diff(y4dot,y4dot);

a81=diff(y4ddot,y1);
a82=diff(y4ddot,y1dot);
a83=diff(y4ddot,y2);
a84=diff(y4ddot,y2dot);
a85=diff(y4ddot,y3);
a86=diff(y4ddot,y3dot);
a87=diff(y4ddot,y4);
a88=diff(y4ddot,y4dot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Jac=[(a11) (a12) (a13) (a14) (a15) (a16) (a17) (a18);
(a21) (a22) (a23) (a24) (a25) (a26) (a27) (a28);
(a31) (a32) (a33) (a34) (a35) (a36) (a37) (a38);
(a41) (a42) (a43) (a44) (a45) (a46) (a47) (a48);
(a51) (a52) (a53) (a54) (a55) (a56) (a57) (a58);
(a61) (a62) (a63) (a64) (a65) (a66) (a67) (a68);
(a71) (a72) (a73) (a74) (a75) (a76) (a77) (a78);
(a81) (a82) (a83) (a84) (a85) (a86) (a87) (a88);]

 

[fresh_42]

I guess the keyword is Wronskian, not Jacobian. Look at
https://www.stewartcalculus.com/dat...ntexts/upfiles/3c3-2ndOrderLinearEqns_Stu.pdf