# Problem with ode15s / ode45

30 views (last 30 days)
Joachim on 29 Jul 2012
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
Patrick is right, specifying time points essentially interpolates the solution obtained by the ode solver, it doesn't necessarily modify the steps taken

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
The method of calling the solver several times is the standard approach to such problems. It is about as "satisfactory" as it gets :) If you want to do just a single call, just write a wrapper function that lets you pass in the pulse times and handles the multiple calls. There shouldn't be too much overhead involved in doing this versus a single call