I need to solve this differential equation using Runge-Kytta 4(5) on Scilab:
The initial conditions are above. The interval and the h-step are:
I don't need to implement Runge-Kutta. I just need to solve this and plot the result on the plane:
I tried to follow these instructions on the official "Scilab Help":
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:
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:
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.
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: