如何在 Scilab 上求解二阶微分方程?
How to solve a second order differential equation on Scilab?
我需要在 Scilab 上使用 Runge-Kytta 4(5) 求解这个微分方程:
初始条件如上。间隔和 h 步长为:
我不需要实施 Runge-Kutta。我只需要解决这个问题并在平面上绘制结果:
我尝试按照官方 "Scilab Help" 上的这些说明进行操作:
建议的代码是:
// 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 的评论,我正在更新此 post。
我需要将二阶 ODE 转换为一阶 ODE 系统,然后我需要编写一个函数来表示这样的系统。
所以,我知道如何用笔和纸来做这件事。因此,使用 z 作为变量:
好的,但是我该如何编写一个普通的脚本?
Xcos 是一次性的。我保留它只是因为我试图模仿官方 Scilab 页面上的示例。
要解决这个问题,你需要使用ode()
,它可以使用很多方法,包括Runge-Kutta。首先,您需要定义一个函数来表示 ODE 系统,您提供的 link 中的 步骤 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();
输出:
我需要在 Scilab 上使用 Runge-Kytta 4(5) 求解这个微分方程:
初始条件如上。间隔和 h 步长为:
我不需要实施 Runge-Kutta。我只需要解决这个问题并在平面上绘制结果:
我尝试按照官方 "Scilab Help" 上的这些说明进行操作:
建议的代码是:
// 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 的评论,我正在更新此 post。 我需要将二阶 ODE 转换为一阶 ODE 系统,然后我需要编写一个函数来表示这样的系统。
所以,我知道如何用笔和纸来做这件事。因此,使用 z 作为变量:
好的,但是我该如何编写一个普通的脚本?
Xcos 是一次性的。我保留它只是因为我试图模仿官方 Scilab 页面上的示例。
要解决这个问题,你需要使用ode()
,它可以使用很多方法,包括Runge-Kutta。首先,您需要定义一个函数来表示 ODE 系统,您提供的 link 中的 步骤 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();
输出: