如何在 Scilab 上求解二阶微分方程？答案

【问题标题】：How to solve a second order differential equation on Scilab?如何在 Scilab 上求解二阶微分方程？
【发布时间】：2019-03-09 04:36:37
【问题描述】：

我需要在 Scilab 上使用 Runge-Kytta 4(5) 求解这个微分方程：

初始条件如上。区间和h-step分别是：

我不需要实现 Runge-Kutta。我只需要解决这个问题并将结果绘制在平面上即可：

我尝试按照官方“Scilab 帮助”上的这些说明进行操作：

https://x-engineer.org/graduate-engineering/programming-languages/scilab/solve-second-order-ordinary-differential-equation-ode-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 cmets，我正在更新这篇文章。我需要将二阶 ODE 转换为一阶 ODE 系统，然后我需要编写一个函数来表示这样的系统。

所以，我知道如何在笔和纸上做到这一点。因此，使用 z 作为变量：

好的，但是我该如何写一个普通的脚本呢？

Xcos 是一次性的。我只是保留它，因为我试图模仿官方 Scilab 页面上的示例。

【问题讨论】：

您展示的代码使用 Xcos 来求解 ODE。那是你要的吗？还是只运行一个常规脚本就足够了？
基本上，你应该做的是将二阶 ODE 转换为一阶 ODE 系统，然后编写一个函数来表示这样的系统。
@luispauloml，谢谢[或者我可以说，obrigado =)]。我更新了我的帖子。我知道如何进行您所说的转换。 XCos 是一次性的，我不需要它。运行常规脚本就足够了。我只需要使用 Runge-Kutta (4,5)
Runge-Kutta-Fehlberg 4(5) 是一种可变步长的方法。是否有关于所需误差容限的信息，最小步长是 0.1，...？
@LutzL 你可以查看ode() help，它有更多关于可用方法中使用的相对和绝对公差的信息。

标签： numerical-methods differential-equations scilab numerical-integration runge-kutta

【解决方案1】：

要解决这个问题，您需要使用ode()，它可以使用多种方法，包括Runge-Kutta。首先，您需要定义一个函数来表示 ODE 系统，您提供的链接中的 Step 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' 以使用任一可用的Runge-Kutta 方法：

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

输出：

【讨论】：