函数lsode可以用于求解形式的ODE
dx -- = f (x, t) dt
使用Hindmarsh的ODE解算器LSODE.
[x, istate, msg] = lsode (fcn, x_0, t) ¶[x, istate, msg] = lsode (fcn, x_0, t, t_crit) ¶常微分方程(ODE)解算器。
要求解的微分方程组为
dx -- = f (x, t) dt
具有
x(t_0) = x_0
解决方案在矩阵中返回x,每行对应于向量的一个元素t。的第一个元素t应该是t_0并且应该对应于系统的初始状态x_0,使得输出的第一行x_0.
第一个参数,fcn,是一个字符串、内联或函数句柄,用于命名函数f调用以计算方程组的右侧向量。函数必须具有以下形式
xdot = f (x, t)
这里的xdot和x是向量和t是标量。
如果fcn是字符串、内联函数或函数句柄的二元素字符串数组或二元素元胞数组,第一个元素命名函数f并且第二元素名称是用于计算的Jacobian的函数f.雅可比函数必须具有以下形式
jac = j (x, t)
这里的jac是偏导数的矩阵
| df_1 df_1 df_1 |
| ---- ---- ... ---- |
| dx_1 dx_2 dx_N |
| |
| df_2 df_2 df_2 |
| ---- ---- ... ---- |
df_i | dx_1 dx_2 dx_N |
jac = ---- = | |
dx_j | . . . . |
| . . . . |
| . . . . |
| |
| df_M df_M df_M |
| ---- ---- ... ---- |
| dx_1 dx_2 dx_N |
第二个参数指定系统的初始状态x_0第三个自变量是向量,t,指定寻求解决方案的时间值。
第四个参数是可选的,可用于指定ODE解算器不应积分的一组时间。这对于解决奇点和导数中存在不连续性的点的困难是有用的。
成功计算后,的值istate将为2(与的Fortran版本一致LSODE).
如果计算不成功,istate会比2更高msg将包含其他信息。
您可以使用该函数lsode_options为设置参数参数lsode.
详见Alan C.Hindmarsh,ODEPACK, A Systematized Collection of ODE Solvers,科学计算,R.S.Stepleman,编辑,(1983)或https://computing.llnl.gov/projects/odepack有关的内部起作用的更多信息lsode.
示例:求解范德波尔方程
fvdp = @(y,t) [y(2); (1 - y(1)^2) * y(2) - y(1)]; t = linspace (0, 20, 100); y = lsode (fvdp, [2; 0], t);
() ¶val = lsode_options (opt) ¶(opt, val) ¶查询或设置函数的参数lsode.
当在没有参数的情况下调用时,将显示所有可用参数的名称及其当前值。
给定一个参数,返回参数的值opt.
当用两个自变量调用时,lsode_options设置参数opt到值val.
参数包括
"absolute tolerance"绝对误差范围。可以是向量,也可以是标量。如果是向量,则它必须与状态向量的维度相匹配。
"relative tolerance"相对误差范围参数。与绝对误差范围不同,此参数只能是标量。
在每个集成步骤中应用的局部错误测试是
abs (local error in x(i)) <= ...
rtol * abs (y(i)) + atol(i)
"integration method"一个字符串,指定用于求解ODE系统的积分方法。有效值为
"adams""non-stiff"没有使用Jacobian(即使它可用)。
"bdf""stiff"使用刚性后向微分公式(BDF)方法。如果没有提供计算雅可比函数,lsode将计算雅可比矩阵的有限差分近似。
"initial step size"在第一个步骤中尝试的步骤大小(默认为自动确定)。
"maximum order"限制求解方法的最大顺序。如果使用Adams方法,此参数必须介于1和12之间。否则,它必须介于1和5之间(包括1和5)。
"maximum step size"设置最大步长将避免经过非常大的区域(未指定默认值)。
"minimum step size"允许的最小绝对步长(默认值为0)。
"step limit"允许的最大步骤数(默认为100000)。
"jacobian type"指定Jacobian类型的字符串,与刚性后向微分公式(BDF)积分方法一起使用。有效值为
"full"默认值。所有的偏导数都是从用户提供的雅可比函数中近似或使用的。
"banded"仅对角线以及参数指定的上下下对角线的数量"lower jacobian subdiagonals"和"upper jacobian
subdiagonals"分别从用户提供的雅可比函数中近似或使用。用户提供的雅可比函数可以将所有其他偏导数设置为任意值。
"diagonal"如果用户提供了雅可比函数,则此设置无效。Jacobian近似为lsode被限制在对角线上,其中每个偏导数都是通过对状态的所有元素一起应用有限变化来计算的;如果真正的雅可比确实总是对角的,这与仅将有限变化应用于状态的相应元素具有相同的效果,但更有效。
"lower jacobian subdiagonals"如果参数,则使用较低的子对角数"jacobian type"设置为"banded"。默认值为零。
"upper jacobian subdiagonals"如果参数,则使用的上部子对角数"jacobian type"设置为"banded"。默认值为零。
以下是使用lsode.给定函数
## oregonator differential equation
function xdot = f (x, t)
xdot = zeros (3,1);
xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) ...
- 8.375e-06*x(1)^2);
xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
xdot(3) = 0.161*(x(1) - x(3));
endfunction
和初始条件x0 = [ 4; 1.1; 4 ],可以使用命令集成一组方程
t = linspace (0, 500, 1000);
y = lsode ("f", x0, t);
如果你这样做,你会查找结果的值在t=0和5,再大约t=305。一组更有效的输出点可能是
t = [0, logspace(-1, log10(303), 150), ...
logspace(log10(304), log10(500), 150)];
上面使用的微分方程的m文件与Octave分布一起包含在examples目录的名称下oregonator.m.
版权所有 © 2024-2025 Octave中文网
ICP备案/许可证号：黑ICP备2024030411号-2