无法使用 Octave 解决简单的 ODE

Question

我是 Octave 的新手，所以在进入更复杂的项目之前，我试图让一些简单的示例起作用。

我正在尝试解析 ODE dy/dx = a*x+b，但没有成功。这是代码：

%Funzione retta y = a*x + b. Ingressi: vettore valori t; coefficienti a,b
clear all;
%Inizializza argomenti
b = 1;
a = 1;
x = ones(1,20);
function y = retta(a, x, b) %Definisce funzione
y = ones(1,20);
y = a .* x .+ b;
endfunction
%Calcola retta
x = [-10:10];
a = 2;
b = 2;
r = retta(a, x, b)
c = b;
p1 = (a/2)*x.^2+b.*x+c  %Sol. analitica di dy/dx = retta %
plot(x, r, x, p1);
% Risolve eq. differenziale dy/dx = retta %
y0 = b; x0 = 0;
p2 = lsode(@retta, y0, x)

输出为：

retta3code
r =

 -18  -16  -14  -12  -10   -8   -6   -4   -2    0    2    4    6    8   10 12   14   16   18   20   22

p1 =

Columns 1 through 18:

82    65    50    37    26    17    10     5     2     1     2     5    10    17    26    37    50    65

Columns 19 through 21:

82   101   122

error: 'b' undefined near line 9 column 16
error: called from:
error:   retta at line 9, column 4
error: lsode: evaluation of user-supplied function failed
error: lsode: inconsistent sizes for state and derivative vectors
error:   /home/fabio/octave_file/retta3code.m at line 21, column 4

因此，函数retta第一次使用正常，但在lsode中使用时失败。为什么会这样？需要更改什么才能使代码正常工作？

Answer 1

下面修改的代码有效，但我更希望能够在函数外定义参数 a 和 b，然后将它们传递给 rdot 作为参数。

x = [-10,10];
a = 1;
b = 0;
c = b;
p1 = (a/2).*(x.^2)+b.*x+c  %Sol. analitica di dy/dx = retta %
function ydot = rdot(ydot, x)
a = 1;
b = 0;
ydot = ones(1,21);
ydot = a.*x .+ b;
endfunction
y0 = p1(1); x0 = 0;
p2 = lsode("rdot", y0, x, x0)'
plot(x, p1, "-k", x, p2, ".r");

Answer 2

不知怎的，你仍然错过了故事的一些重要部分。要求解 ODE y'=f(y,x)，您需要定义一个函数

function ydot = f(y,x)

其中 ydot 与 y 具有相同的维度，两者都必须是向量，即使它们是维度 1。x 是标量。由于某些传统原因，lsode（在多个求解器包中使用的 FORTRAN 代码）更喜欢使用较少的顺序 (y,x)，在大多数教科书和其他求解器中，您会找到顺序 (x,y).

然后为了在样本点 xlist 上获得解决方案样本 ylist 你调用

ylist = lsode("f", y0, xlist)

其中 xlist(1) 是初始时间。

f 的内部结构与示例列表 list 及其大小无关。这是一个单独的问题，您可以使用多重评估来计算

之类的确切解决方案

yexact = solexact(xlist)

给pass parameters, use anonymous functions，喜欢

function ydot = f(y,x,a,b)
    ydot = [ a*x+b ]
end

a_val = ...
b_val = ...
lsode(@(y,x) f(y,x,a_val, b_val), y0, xlist)

无法使用 Octave 解决简单的 ODE

Failing to solve a simple ODE with Octave

function

octave

ode