How to solve a second order ordinary differential equation (ODE) in Scilab

In the tutorial How to solve an ordinary differential equation (ODE) in Scilab we can see how a first order ordinary differential equation is solved (numerically) in Scilab. In this tutorial we are going to solve a second order ordinary differential equation using the embedded Scilab function ode().

As example we are going to use a nonlinear second order ordinary differential equation:

\[\frac{d^2 y}{dt^2} = \frac{1}{t+1} + sint(t)\sqrt{t}\]

with the initial conditions:

\[ \begin{equation*} \begin{split}
y(0) &= 0\\
\frac{dy}{dt}(0) &= -2
\end{split} \end{equation*} \]

The first step is to write the second order differential equation as a system of two first order differential equations:

\[ \begin{equation*} \begin{split}
\frac{dy}{dt} &= z\\
\frac{dz}{dt} &= \frac{1}{t+1} + sint(t)\sqrt{t}
\end{split} \end{equation*} \]

This way we have set our differential equation in the form:

\[\frac{dx}{dt}=f(t,x)\]

where both dx and x are matrices of size 2 x 1.

\[ \begin{equation*} \begin{split}
x &= \left [ z; \frac{1}{t+1} + sint(t)\sqrt{t} \right ]\\
dx &= \left [ \frac{dy}{dt}; \frac{dz}{dt} \right ]
\end{split} \end{equation*} \]

To solve this differential equation in Scilab, first we need to define our differential equation as a separate function. Scilab allows to define a custom function is an *.sce file, together with other instructions. For this example, all of the Scilab instruction will need to be included in the same *.sce file.

Step 1. Define the differential equation as a custom function

function dx = f(t, x)
 dx(1) = x(2);
 dx(2) = 1/(t+1) + sin(t)*sqrt(t);
endfunction

Step 2. Define the integration time t, initial time t₀ and initial values

t = 0:0.01:5*%pi;
t0 = min(t);
y0 = [0; -2];

Step 3. Call the ode() function with the above parameters

y = ode(y0, t0, t, f);

Step 4. Plot the result. Notice that the numerical solution y contains the second and first integration results, y and y’.

plot(t,y(1,:),'LineWidth',2)
plot(t,y(2,:),'r','LineWidth',2)
xgrid();
xlabel('$t \quad [s]$','FontSize',3)
ylabel('$f(t,x)$','FontSize',3)
title(['Integration of ' '$\frac{d^2 x}{dt^2} = \frac{1}{t+1} + sint(t)\sqrt{t}$'],'FontSize',3)
legend(['$\Large{x}$' '$\Large{dx/dt}$'],2)

By running all of the above Scilab instruction in an script file (*.sce), we get the following graphical window:

To check that the solution of our integration is correct, we are going the model the equation in Xcos and run the simulation for 15.71 seconds (5π).

The Xcos block diagram model of the second order ordinary differential equation is integrated using the Runge-Kutta 4(5) numerical solver.

After running the simulation, Xcos will output the following graphical window (the grid has been added afterwards):

As you can see, both methods give the same results. This is a confirmation that the system of first order ODE were derived correctly and the equations were correctly integrated.

For any questions, observations and queries regarding this article, use the comment form below.

Don’t forget to Like, Share and Subscribe!