[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Using Hindmarsh’s ODE solver LSODE in OCTAVE

From: Sebastian Schöps
Subject: Re: Using Hindmarsh’s ODE solver LSODE in OCTAVE
Date: Mon, 28 Dec 2015 04:05:22 -0800 (PST)

Omer Tzuk wrote
> I have a question on using LSODE package in OCTAVE, I have put it on
> stackoverflow -
> I will be grateful if you can refer to my question,

Crossposting is not very welcome
( and it's quite an effort that
someone has to make to help you: I have to go from here to stackoverflow and
then to github and finally, your code is not commented at all. What kind of
system do you want to solve? I suggest to write it down first.

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


> T = 0:1e-2:1; % time vector
> K = sprandsym(32,1)+eye(32,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.


View this message in context:
Sent from the Octave - General mailing list archive at

reply via email to

[Prev in Thread] Current Thread [Next in Thread]