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