Solving a system of second order tightly coupled nonlinear ODE with six initial conditions in Matlab

Question

I am solving a problem from fluid dynamics; in particular tightly coupled nonlinear ordinary differential equations. The following is a scaled-down version of my actual problem.

I have solved system of coupled odes many times in the past but this case is different since double derivatives of one variable depends on the double derivative of another variable.

How do I implement it in $ode45$? I need $3$ x $2 = 6$ plots of $x$, $x-dot$ and $x-ddot$ versus time for $t$, $0$ to $2$. All required initial conditions have zero values.

https://postimg.cc/0zf8XV9x

Any help will be highly appreciated.
If the ODEs had no last terms of double derivatives, I would have done it myself.

Answer 1

Your system linked to in the question, \begin{align} \ddot x_1&=t^2x_1\sin(4\pi t)-x_2^2-\dot x_1x_2t+\dot x_2t+\ddot x_2\\ \ddot x_2&=t^2x_2\sin(4\pi t)-x_1^2+\dot x_1x_2t+\dot x_1t-\ddot x_1 \end{align} is just a linear system for the second derivatives, with constant system matrix and the right side depending on the values and first derivatives. In the example $$ \pmatrix{1&-1\\1&1}\pmatrix{\ddot x_1\\\ddot x_2}=\pmatrix{t^2x_1\sin(4\pi t)-x_2^2-\dot x_1x_2t+\dot x_2t\\t^2x_2\sin(4\pi t)-x_1^2+\dot x_1x_2t+\dot x_1t}. $$ Construct the matrix, or get it as global constant, construct the right side, and solve the system with one of the provided linear solvers.

def odefun(t, u):
    x1, x2, dotx1, dotx2 = u
    A = np.array([[1,-1],[1,1]]);
    s = np.sin(4*np.pi*t)*t**2
    b = np.array([ 
         s*x_1-x_2**2-dotx1*x2*t+dotx2*t, 
         s*x_2-x_2**2+dotx2*x1*t+dotx1*t
    ])
    dot2x = np.linalg.solve(A,b);
    return [ dotx1, dotx2, dot2x[0], dot2x[1] ] 

With the matlab solvers, you can this even leave un-solved and provide the system matrix of the full sytem $$ \pmatrix{1&0&0&0\\0&1&0&0\\0&0&1&1\\0&0&-1&1} \pmatrix{\dot x_1\\\dot x_2\\\ddot x_1\\\ddot x_2} = \pmatrix{\dot x_1\\\dot x_2\\b_1(t,x_1,x_2,\dot x_1,\dot x_2)\\b_2(t,x_1,x_2,\dot x_1,\dot x_2)}. $$ as the weight/mass matrix parameter.

---

The modified system that you produced in the comments, \begin{align} \ddot x_1&=tx_2-x_1\dot x_2-\ddot x_2\\ \ddot x_2&=t^2x_1+x_2\dot x_1-\ddot x_1 \end{align} looks more simple but is in fact more complicated as it is singular in the higher order derivatives. Indeed, one can replace the second equation be the difference of both to get \begin{align} \ddot x_1+\ddot x_2&=tx_2-x_1\dot x_2\\ 0&=t^2x_1-tx_2+x_2\dot x_1+x_1\dot x_2 \end{align} This can be solved into an explicit first order system by introducing the state vector $(x_1,x_2,u)$ with the intended relation $u=\dot x_1+\dot x_2$ so that \begin{align} \dot x_1&=\frac{t^2x_1-tx_2+ux_1}{x_1-x_2},\\ \dot x_2&=\frac{t^2x_1-tx_2+ux_2}{x_2-x_1},\\ \dot u&=tx_2-x_1\frac{t^2x_1-tx_2+ux_2}{x_2-x_1} =\frac{tx_2^2-t^2x_1^2-ux_1x_2}{x_2-x_1}. \end{align} An actual implementation could re-use intermediary expressions like the value of $\dot x_2$.