help-octave
[Top][All Lists]
Advanced

[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 -
> http://stackoverflow.com/q/34453043/2051572
> 
> I will be grateful if you can refer to my question,

Crossposting is not very welcome
(https://en.wikipedia.org/wiki/Crossposting) 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);

Then

> 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.

Bye
Sebastian



--
View this message in context: 
http://octave.1599824.n4.nabble.com/Using-Hindmarsh-s-ODE-solver-LSODE-in-OCTAVE-tp4674210p4674211.html
Sent from the Octave - General mailing list archive at Nabble.com.



reply via email to

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