# Numerically solving a system of linear 2nd order differential equations

Background

I am working on a project for differential equations, and I need to demonstrate to the professor that I am capable of 1) learning new techniques on my own and 2) implementing those techniques in a programming environment (MATLAB, in this case).

My question concerns how to solve a 2nd order system of differential equations using numerical methods. If someone wants to provide a full answer or a sketch of the solution, I would be very happy! Otherwise, just pointing me in the right direction, perhaps to a particular method, website, or book, would be helpful.

Note: While my system of equations can be solved symbolically, the project explicitly calls for the use of numerical methods. I CANNOT use the numerical functions already built into MATLAB. I have to build my own using the appropriate numerical methods.

The System of Equations

I have the following linear 2nd order system of differential equations:

$x''_1 = 4x_1 + 2x_2 + 2cos(3t) \\x''_2 = -4x_2 + 2x_1 + 2x_3 \\x''_3 = -2x_3 + 2x_2$

which has initial conditions $x_1(0) = 0, x_2(0) = 0, x_3(0) = 0, x'_1(0) = 0, x'_2(0) = 0, x'_3(0) = 0$.

The system can be rewritten as a matrix, as seen below.

$\left[ \begin{matrix} x_1'' \\ x_2'' \\ x_3'' \\ \end{matrix} \right]=$ $\left[ \begin{matrix} -4&2&0 \\ 2&-4&2 \\ 0&2&-2 \\ \end{matrix} \right]$ $\left[ \begin{matrix} x_1 \\ x_2 \\ x_3 \\ \end{matrix} \right]$ $+ \left[ \begin{matrix} 2cos(3t) \\ 0 \\ 0 \\ \end{matrix} \right]$

In addition, we can transform this system of 2nd order equations into a system of 1st order equations. Transforming the original three equations yields

$a_1 = x_1 \\ a_2 = x_1' \\ b_1 = x_2 \\ b_2 = x_2' \\ c_1 = x_3 \\ c_2 = x_3'$

$a_1' = a_2 \\ a_2' = -4a_1 + 2b_1 + 2cos(3t) \\ b_1' = b_2 \\ b_2' = -4b_1 + 2a_1 + 2c_1 \\ c'_1 = c_2 \\ c'_2 = -2c_1 + 2b_1$

Now what?

At this point, I'm not sure what to do. I have used Newton's method to approximate simple functions before, but not 2nd order DE's and certainly not a system of DE's.

What methods would be appropriate here? How should I proceed? Any help would be greatly appreciated.

• Do these equations come with initial conditions or with boundary conditions? The problem is very different in these two cases. My guess from the context would be that you have initial conditions. In this case the most straightforward thing to do are explicit methods. The simplest of these is Euler's method, but this method has severe problems (very small step sizes are required for small errors at each step, and sometimes even smaller step sizes are required for stability). Two others are the explicit Runge midpoint rule and the explicit Runge trapezoidal rule. – Ian Mar 17 '16 at 19:07
• An explicit method which is popular because it nicely balances out various things is classical 4th order Runge Kutta. For a low dimensional problem like this one, most numerical analysts would probably start with the Dormand-Price method and see if it works, in part because it is implemented in Matlab's ode45 function. But you could implement it yourself without too much pain. – Ian Mar 17 '16 at 19:07
• Implicit methods avoid some of the stability problems associated with explicit methods, but they come at the cost of having to solve a system of algebraic equations at each time step. The simplest one is backward Euler, but there are higher order ones as well. – Ian Mar 17 '16 at 19:09
• @Ian They do come with initial conditions. I added them to my original question (above). Is there an easy way to explain the jump from single equations to a system of equations? How does the implementation of these methods change? – EternusVia Mar 17 '16 at 19:12
• Once you've recast into a first order system, everything about the methods is the same, you are just adding vectors instead of adding scalars. – Ian Mar 17 '16 at 19:12