Search code examples
numerical-methodsdifferential-equationsscilabnumerical-integrationrunge-kutta

How to solve a second order differential equation on Scilab?


I need to solve this differential equation using Runge-Kytta 4(5) on Scilab:

enter image description here

The initial conditions are above. The interval and the h-step are:

enter image description here

enter image description here

I don't need to implement Runge-Kutta. I just need to solve this and plot the result on the plane:

enter image description here

I tried to follow these instructions on the official "Scilab Help":

https://x-engineer.org/graduate-engineering/programming-languages/scilab/solve-second-order-ordinary-differential-equation-ode-scilab/

The suggested code is:

// Import the diagram and set the ending time
loadScicos();
loadXcosLibs();
importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.zcos");
scs_m.props.tf = 5000;

// Select the solver Runge-Kutta and set the precision
scs_m.props.tol(6) = 6;
scs_m.props.tol(7) = 10^-2;

// Start the timer, launch the simulation and display time
tic();
try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
t = toc();
disp(t, "Time for Runge-Kutta:");

However, it is not clear for me how I can change this for the specific differential equation that I showed above. I have a very basic knowledge of Scilab.

The final plot should be something like the picture bellow, an ellipse:

enter image description here

Just to provide some mathematical context, this is the differential equation that describes the pendulum problem.

Could someone help me, please?

=========

UPDATE

Based on @luizpauloml comments, I am updating this post. I need to convert the second-order ODE into a system of first-order ODEs and then I need to write a function to represent such system.

So, I know how to do this on pen and paper. Hence, using z as a variable:

enter image description here

OK, but how do I write a normal script?

The Xcos is quite disposable. I only kept it because I was trying to mimic the example on the official Scilab page.


Solution

  • To solve this, you need to use ode(), which can employ many methods, Runge-Kutta included. First, you need to define a function to represent the system of ODEs, and Step 1 in the link you provided shows you what to do:

    function z = f(t,y)
        //f(t,z) represents the sysmte of ODEs:
        //    -the first argument should always be the independe variable
        //    -the second argument should always be the dependent variables
        //    -it may have more than two arguments
        //    -y is a vector 2x1: y(1) = theta, y(2) = theta'
        //    -z is a vector 2x1: z(1) = z    , z(2) = z'
    
        z(1) = y(2)         //first equation:  z  = theta'
        z(2) = 10*sin(y(1)) //second equation: z' = 10*sin(theta)
    endfunction
    

    Notice that even if t (the independent variable) does not explicitly appear in your system of ODEs, it still needs to be an argument of f(). Now you just use ode(), setting the flag 'rk' or 'rkf' to use either one of the available Runge-Kutta methods:

    ts = linspace(0,3,200);
    theta0  = %pi/4;
    dtheta0 = 0;
    y0 = [theta0; dtheta0];
    t0 = 0;
    thetas = ode('rk',y0, t0, ts, f); //the output have the same order
                                      //as the argument `y` of f()
    
    scf(1); clf();
    plot2d(thetas(2,:),thetas(1,:),-5);
    xtitle('Phase portrait', 'theta''(t)','theta(t)');
    xgrid();
    

    The output:

    Phase portrait of a gravity pendulum