MATLAB: Problem with ode15s / ode45

differential equations

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

Best Answer

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