30 views (last 30 days)

Show older comments

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.

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.

Richard Brown
on 10 Aug 2012

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

Start Hunting!