在 OCTAVE 中使用 Hindmarsh 的 ODE 求解器 LSODE

Question

正在学习OCTAVE，正在尝试使用LSODE ODE solver集成一个版本FitzHugh–Nagumo model。我的尝试是这样的：

time = linspace(0,200,1000);
u0 = rand(32,32);
v0 = rand(32,32);

vec_u0 = reshape(u0,[1,size(u0)(1)*size(u0)(2)]);
vec_v0 = reshape(v0,[1,size(v0)(1)*size(v0)(2)]);
vec_FHN0 = horzcat(vec_u0,vec_v0);

FHN = lsode("FHN_eq_vec", vec_FHN0, time);
FHN(end)

我定义的所有函数都在我在 GitHub 中设置的存储库中 - link. I have created a function that transform the two 2D fields of the FHN model into a row vector (as I understand from the examples 这里 LSODE 集成器使用行向量作为输入）。我收到此错误消息：

>> fhn_integrate_lsode
warning: non-integer range used as index
warning: called from
    FHN_eq_vec at line 3 column 7
    fhn_integrate_lsode at line 9 column 5
error: reshape: can't reshape 0x1 array to 1x1 array
error: called from
    FHN_eq_vec at line 4 column 3
    fhn_integrate_lsode at line 9 column 5
error: lsode: evaluation of user-supplied function failed
error: called from
    fhn_integrate_lsode at line 9 column 5
error: lsode: inconsistent sizes for state and derivative vectors
error: called from
    fhn_integrate_lsode at line 9 column 5
>>

有人知道可能是什么问题吗？

Answer 1

已在 http://octave.1599824.n4.nabble.com/Using-Hindmarsh-s-ODE-solver-LSODE-in-OCTAVE-td4674210.html

回答此问题

However, looking quickly at your code, the system that you want to solve is probably an ordinary differential equation stemming from a pde space discretization, i.e.,

$dx(t)/dt = f(x,t) := -K x(t) + r(t)$

with K being a square matrix (Laplacian?!) and f a time-dependent function of matching dimension. I expect that your system is stiff (due to the negative Laplacian on the right-hand side) and that you are happy with errors in the order of 10^(-4). Thus you should adapt the options of lsode:
lsode_options("integration method","stiff");
lsode_options("absolute tolerance",1e-4);
lsode_options("relative tolerance",1e-4);
Then
T = 0:1e-2:1; % time vector
K = sprandsym(32,1)+eye(32); % symmetric stiffness matrix
x0 = rand(32,1); % initial values
r = @(t) rand(32,1)*sin(t); % excitation vector
f = @(x,t) (-K*x+r(t)); % right-hand-side function
x=lsode (f, x0, T); % get solution from lsode
You should exploit any knowledge on the Jacobian df/dx since this will speed up computations. It's trivial in the case of linear ODE:
f = {@(x,t) -K*x+r(t), @(x,t) -K}; % right-hand-side function.
On the other hand, if the system has an additional mass matrix

$M dx(t)/dt = -K x(t) + r(t)$

things might become more complicated. You probably want to use another time stepper. Only if M has full rank then you could do
f = @(x,t) ( M\(-K*x+r(t)) ); % right-hand-side function
which is typically not very efficient.

Bye Sebastian

在 OCTAVE 中使用 Hindmarsh 的 ODE 求解器 LSODE

Using Hindmarsh’s ODE solver LSODE in OCTAVE

matlab

octave

numerical-integration