79 views (last 30 days)

Hi together, I'm very new to solving differential equations in Matlab... Here I face the following problem. I have a (simplified) model with two states, where an external series of pulses drives the population from one state to the other. (the output figure should explain the situation).

My problem is, that the code works fine when I set the 'MaxStep' property to 0.1. If I do not set the property or use higher values, the ODEsolver misses all the drivig pulses. This may be due to the fact, that the whole calculating time-interval is 8000 time units and the pulses are sparse and short (like a 0.1-long pulse every 200 time units).

This situation is kind of unsatisfactory. Is there a certain rule how to set the MaxStep property? Is there any way to tell the solver in which timeintervals the pulses are (where the stepsize would then have to be reduced)? Or to "presample" the driving function?

Any suggestions? Thank you in advance.

Joachim

function run

%parameters for the driving function

exc=1; %pulse amplitude

t_p=10; %time of first pulse (ns)

sigma_p=0.02; %pulse width in ns (sigma)

exc_burst=[0 1000]; %duration of excitation burst in ns: [start stop]

repetition_cycle=200; %pulse-to-pulse period of the pulse trains in ns.

%driving Function (series of Gaussian pulses).

function erg=I_drv(t)

erg=exc/sqrt(2)/pi/sigma_p*exp(-(((t-t_p)-(floor(((t-t_p)+repetition_cycle/2)/repetition_cycle)*repetition_cycle)).^2)./2/(sigma_p.^2));

end

%ODE set

function [Ndot] = myode(t,N)

% S0

Ndot(1) = (-I_drv(t).^2)*0.0001*N(1);

% S1

Ndot(2) = (+I_drv(t).^2)*0.0001*N(1);

Ndot = Ndot'; % This makes Ydot into a column vector

end

%solve ODEs

A=0; %start time in ns

B=20000; %stop time

options = odeset('RelTol', 2.22045e-14,'AbsTol',[1e-16 1e-16],'MaxStep',0.1);

% options = odeset('RelTol', 2.22045e-14,'AbsTol',[1e-16 1e-16]);

[T,Y]=ode15s(@myode,[A B],[1 0],options);

%Plot ODEs

f_handle=figure(1);

set(f_handle,'position',[600 50 670 870]);

clf

subplot(2,1,1)

title 'driving pulses'

cla

hold on

plot(A:0.001:B,I_drv(A:0.001:B),'red');

subplot(2,1,2)

title 'state populations'

cla

plot(T,Y)

legend('N1', 'N2')

ylim([0 1])

end

Star Strider
on 29 Jul 2012

Rather than specifying a two-element vector for ‘tspan’, I suggest you give it a vector of the times at which you want it to integrate your equations. From the documentation for ‘ode23’ and the other ODE functions (I put the relevant sections in bold):

--------------------------------------------

tspan —— A vector specifying the interval of integration, [t0,tf]. The solver imposes the initial conditions at tspan(1), and integrates from tspan(1) to tspan(end). To obtain solutions at specific times (all increasing or all decreasing), use tspan = [t0,t1,...,tf].

For tspan vectors with two elements [t0 tf], the solver returns the solution evaluated at every integration step. For tspan vectors with more than two elements, the solver returns solutions evaluated at the given time points. The time values must be in order, either increasing or decreasing.

--------------------------------------------

I do this frequently, because it also shortens the time it takes the ODE functions to integrate my ODEs. The ‘tspan’ vector elements do not have to be uniformly-spaced, simply non-decreasing.

Patrick
on 1 Aug 2012

Actually, If I am correct, specifying tspan does not change the result (see next paragraphs in documentation):

--------------------------------------------

Specifying tspan with more than two elements does not affect the internal time steps that the solver uses to traverse the interval from tspan(1) to tspan(end). All solvers in the ODE suite obtain output values by means of continuous extensions of the basic formulas. Although a solver does not necessarily step precisely to a time point specified in tspan, the solutions produced at the specified time points are of the same order of accuracy as the solutions computed at the internal time points.

Specifying tspan with more than two elements has little effect on the efficiency of computation, but for large systems, affects memory management.

--------------------------------------------

My guess is that the solver does not naturally "see" the pulses. Therefore, if the time steps do not happen to hit them, they are missed. Question remains: how to tell the solver to hit specific time steps? An unsatisfactory answer is to call it several times from a loop ; in that case, it can be called "when there is a pulse" and "when there is no pulse" but this is probably time consuming.

Richard Brown
on 10 Aug 2012

Jan
on 30 Jul 2012

Edited: Jan
on 30 Jul 2012

'RelTol', 2.22045e-14, 'AbsTol', [1e-16 1e-16]

This is ridiculously small. It is very unlikely that you get any meaningful data with such small tolerances: Because the 64-bit doubles allow a relative accuracy of about 1e-16, the relative tolerance of 2e-14 means (roughly spoken), that that your calculations use 2 significant digits only.

It is recommended to scale the problem, such that all values are near to 1. Then use large tolerances at first, e.g. 1e-2 relative and 1e-4 absolute. Afterwards you can reduces the tolerances and observe the effects to the results.

Patrick
on 3 Aug 2012

Have you tried to use events or to call the solver several times as I suggested in my comment above?

Richard Brown
on 10 Aug 2012

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.