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`

.

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 _{0}* 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!

## Hynek Pokorný

Hello,

I do not understand the function deffinition. Why do you put down dx(1) = x(2)? The corresponding value to “dx(1)” is “z”, which is on the first place in the vector, therefore I would expect “x(1)”.

Can somebody explain me? Thanks, Hynek