Failing to solve a simple ODE with Octave

Question

I am new to Octave, so I am trying to make some simple examples work before moving onto more complex projects.

I am trying to resolve the ODE dy/dx = a*x+b, but without success. Here is the code:

%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)

And the output is:

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

So, the function retta works properly the first time, but it fails when used in lsode. Why does that happen? What needs to be changed to make the code work?

Answer 1

Somehow you still miss some important parts of the story. To solve an ODE y'=f(y,x) you need to define a function

function ydot = f(y,x)

where ydot has the same dimensions as y, both have to be vectors, even f they are of dimension 1. x is a scalar. For some traditional reason, lsode (a FORTRAN code used in multiple solver packages) prefers the less used order (y,x), in most text books and other solvers you find the order (x,y).

Then to get solution samples ylist over sample points xlist you call

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

where xlist(1) is the initial time.

The internals of f are independent of the sample list list and what size it has. It is a separate issue that you can use multi-evaluation to compute the exact solution with something like

yexact = solexact(xlist)

---

To pass parameters, use anonymous functions, like in

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)