Как решить дифференциальное уравнение второго порядка на Scilab?

Question

Мне нужно решить это дифференциальное уравнение, используя Runge-Kytta 4 (5) на Scilab:

Начальные условия выше.Интервал и h-шаг:

Мне не нужно реализовывать Рунге-Кутта.Мне просто нужно решить эту проблему и отобразить результат на плоскости:

Я пытался следовать этим инструкциям на официальной «Справке Scilab»:

https://x -engineer.org / аспирантура / языки программирования / scilab / решение второго порядка обыкновенное-дифференциальное уравнение-од-scilab /

Предлагаемый код:

// 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:");

Однако мне не ясно, как я могу изменить это для конкретного дифференциального уравнения, которое я показал выше.У меня есть базовые знания о Scilab.

Финальный сюжет должен выглядеть примерно так, как показано на рисунке ниже, эллипс:

Просто чтобы дать некоторый математический контекст, это дифференциальное уравнение, которое описывает проблему маятника.

Может кто-нибудь помочь мне, пожалуйста?

=========

ОБНОВЛЕНИЕ

Основываясь на комментариях @luizpauloml, я обновляю этот пост.Мне нужно преобразовать ODE второго порядка в систему ODE первого порядка, а затем мне нужно написать функцию для представления такой системы.

Итак, я знаю, как это сделать на ручке и бумаге.Следовательно, используя z в качестве переменной:

ОК, но как мне написать обычный скрипт?

Xcos вполнеодноразовые.Я сохранил его только потому, что пытался подражать примеру на официальной странице Scilab.

Answer 1

Чтобы решить эту проблему, вам нужно использовать ode(), который может использовать много методов, включая Рунге-Кутту.Во-первых, вам нужно определить функцию для представления системы ODE, и Шаг 1 в предоставленной вами ссылке покажет вам, что делать:

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

Обратите внимание, что даже если t (независимая переменная) явно не появляется в вашей системе ODE, она все еще должна быть аргументом f().Теперь вы просто используете ode(), устанавливая флаг 'rk' или 'rkf', чтобы использовать любой из доступных методов Рунге-Кутты:

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();

Вывод: