Solving a system of differential equations and plotting solutions with Scilab

Question

I'm trying to solve this system of differential equations with Scilab :

It describes the movement of a marble sliding on a plane. I would like my Scilab programme to plot the trajectory of the marble, i.e. rcos(theta) on the X-axis and rsin(theta) on the Y-axis.

I'm completely new to Scilab, I've never had any course on how to use it yet. I wrote the following code based on some examples given by my teacher. So feel free to explain like I'm 5.

Here's my code :

// constants definition

m1 = 1 ;
m2 = 1 ;
v0 = 1 ; 
l = 1 ; 
g = 9.81 ;

function dxdt = f(t,x)
    // returns [r ddot ; theta dot]
    r = x(1) ;
    theta = x(2) ;
    
    r_ddot = m1*(l*v0)^2/(r^3*(m1+m2)) - m2*g ;// r ddot
    theta_dot = l*v0/r^2 ;// θ dot
    
    dxdt = [r_ddot; theta_dot] ;
endfunction

// initial conditions definition

r0 = 0 ;
theta0 = 0 ;

x0 = [r0,theta0] ;

// simulation range definition

t0 = 0 ; 
tf = 20 ;
Nb_points = 500 ;

time=[t0:(Nb_points-1)/(tf-t0):tf];

// system resoution

Evo_x=ode(x0,t0,time,f);

// graph plotting

r = Evo_x(:, 1) ;
theta = Evo_x(:, 2) ;

xdata = r.*cos(theta) ;
ydata = r.*sin(theta) ;

plot(xdata, ydata);
xlabel("x") ;
ylabel("y") ;
title ("Trajectory of the M1 marble")

It compiles succesfully, but the trajectory doesn't appear on the resulting graph. empty graph

Can you help me spot the problems ? Thank you very much in advance.

Answer 1

You had some problems in your program. First, r_ddot in the right-hand-side was not correct w.r.t your equations (I fixed it below). Moreover you forgot one equation for r as the equation is second order. The state of your first order equivalent equation is [r,r_dot,theta]. Last problem, you should not initialize r to 0 (r is in the denominator) and your time vector was not correctly defined.

// constants definition

m1 = 1 ;
m2 = 1 ;
v0 = 1 ; 
l = 1 ; 
g = 9.81 ;

function dxdt = f(t,x)
    // returns [r dot; r ddot ; theta dot]
    r = x(1) ;
    r_dot = x(2);
    theta = x(3) ;
    
    r_ddot = (m1*(l*v0)^2/r^3 - m2*g)/(m1+m2) ;// r ddot
    theta_dot = l*v0/r^2 ;// θ dot
    
    dxdt = [r_dot; r_ddot; theta_dot] ;
endfunction

// initial conditions definition

r0 = 1 ;
rd0 = 0;
theta0 = 0 ;

x0 = [r0;rd0;theta0] ;

// simulation range definition

t0 = 0 ; 
tf = 20 ;
Nb_points = 500 ;

time=linspace(t0,tf,Nb_points);

// system resolution

Evo_x=ode(x0,t0,time,f);

// graph plotting

r = Evo_x(1,:) ;
theta = Evo_x(3,:) ;

xdata = r.*cos(theta) ;
ydata = r.*sin(theta) ;

plot(xdata, ydata);
xlabel("x") ;
ylabel("y") ;
title ("Trajectory of the M1 marble")