How to solve an ordinary differential equation (ODE) in Scilab

Scilab comes with an embedded function for solving ordinary differential equations (ODE). For a better understanding of the syntax we are going to solve an ODE analytically. For example, let’s have a look at the following ordinary differential equation:

\[y\prime=\frac{x+1}{y},\quad y(0)=0.1\]

where:

y – is the dependent variable (the equation contains the derivative of y)
x – is the independent variable (the derivative is with respect to x)
y(0)=0.1 – is the initial condition, at x = 0 , y = 0.1

We will solve this differential equation analytically. First we’ll write it again, this time using the Leibniz notation (dy/dx).

\[\frac{dy}{dx}=\frac{x+1}{y},\quad y(0)=0.1\]

The equation is with separable variables and it’s solved as:

\[ \begin{equation*} \begin{split}
\frac{dy}{dx} &= \frac{x+1}{y}\\
ydy &= (x+1)dx\\
\int ydy &= \int (x+1)dx + C\\
\frac{y^2}{2} &= \frac{x^2}{2}+x+C\\
y^2 &= x^2+2x+C\\
y &= \sqrt{x^2+2x+C}
\end{split} \end{equation*} \]

The general solution of our differential equation is:

\[y(x)=\sqrt{x^2+2x+C}\]

To find the value of the constant C, we need to use the initial condition:

\[ \begin{equation*} \begin{split}
y(0) &= 0.1\\
\sqrt{0^2+2 \cdot 0+C} &= 0.1\\
C &= 0.1^2\\
C &= 0.01
\end{split} \end{equation*} \]

Now we have our particular solution of the differential equation:

\[y(x)=\sqrt{x^2+2x+0.01}\]

We will use this solution to compare against the result of the numerical integration.

In order to solve in Scilab an ordinary differential equation, we can use the embedded function ode(). This simplest syntax of the function is:

Scilab syntax for ode() function

Image: Scilab syntax for ode() function

where:

y – is the return (dependent) variable, the solution of the differential equation; it can be a vector or a matrix, depending on the number of differential equations
y0 – is the initial condition of the differential equation; can be a real vector or matrix
x0 –  is the initial value of the independent variable; is a real scalar
x – a real vector, the values of the independent variable for which the solution is calculated
f – is a function, external, string or list, representing the right hand side of the differential equation

The ode() function invokes a numerical method, which solves the differential equation numerically. By default lsoda solver of package ODEPACK is called. It automatically selects between nonstiff predictor-corrector Adams method and stiff Backward Differentiation Formula (BDF) method. It uses nonstiff method initially and dynamically monitors data in order to decide which method to use.

To find the numeric solution, first, we need to define our differential equation. We’ll do this by using the Scilab function deff():

deff('yprim=f(x,y)','yprim=(x+1)/y');

Second, we will define the values of x for which we want to compute the solution of the differential equation. We will choose x between 0 and 1 with an increment of 0.001.

x0=0; xinc=0.001; xf=1; x=x0:xinc:xf;

Third, we write the initial value of the solution and we call the ode() function with the appropriate parameters.

y0=0.1;
ydiff=ode(y0,x0,x,f);

In order to proceed faster with the tutorial, you can copy and paste into an *.sce file, the Scilab instructions below. The first part is the evaluation of the analytic solution, the second part is the computation of the numeric solution, using the ode() function.

// Define x
x0=0; xinc=0.001; xf=1; x=x0:xinc:xf;
// Calculate analytic solution
y=sqrt(x.^2+2*x+0.01);
// Plot analytic solution
subplot(2,1,1), plot(x,y), xgrid
ylabel('y(x)','fontsize',2)
title('Analytic solution','fontsize',2)

// Define differential equation
deff('yprim=f(x,y)','yprim=(x+1)/y');
// Solve differential equation
y0=0.1;
ydiff=ode(y0,x0,x,f);
// Plot numeric solution
subplot(2,1,2), plot(x,ydiff,'r'), xgrid
title('Numeric solution','fontsize',2)
ylabel('y(x)','fontsize',2)
xlabel('x','fontsize',2)

The result of the analytic solution is the variable y. The result of the numeric solution is the variable ydiff. The two variables are plotted in a graphical window, on different axes. You can plot both y and ydiff on the same axis but, for this particular example, the values overlapped (which means that the accuracy of the numerical method is very good).

Analytic and numeric solution for an ODE

Image: Analytic and numeric solution for an ODE

The equation above was a linear ordinary differential equation. Let’s use the ode() function to solve a nonlinear ODE.

\[y\prime=y^2-\sqrt{t},\quad y(0)=0\]

Notice that the independent variable for this differential equation is the time t. The solution as well as the graphical representation are summarized in the Scilab instructions below:

t0=0; tf=4; dt=0.001; t=t0:dt:tf;
deff('yprim=f(t,y)','yprim=y^2-sqrt(t)');
y0=0;
y=ode(y0,t0,t,f);
plot(t,y), xgrid
xlabel('$x$','fontsize',4)
ylabel('$y(x)$','fontsize',4)
title('$y\prime=y^2-\sqrt{t},\quad y(0)=0$','fontsize',4)

By running the above Scilab instructions, we get the graphical representation of the numeric solution y:

Numeric solution for a nonlinear differential equation

Image: Numeric solution for a nonlinear differential equation

The function ode() is very flexible and allows the user to choose the numerical integration method and also the relative and absolute tolerances (error). How to use these settings, we’ll discuss in detail in a different tutorial.

For any questions or observations regarding this tutorial please use the comment form below.

Don’t forget to Like, Share and Subscribe!

One Response

  1. Dr K V B RAJAKUMAR

Leave a Reply

Ad Blocker Detected

Dear user, Our website provides free and high quality content by displaying ads to our visitors. Please support us by disabling your Ad blocker for our site. Thank you!

Refresh