Problem with ode15s / ode45

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

Answer 1

[Star Strider]

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.

[Joachim]

Hi there,

thank you for your fast reply.

This sounds like exactly what I need... However, I cannot get it working properly. Even when I use tspan=[A:0.001:B] Matlab misses all the pulses and does not output the correct solution.

Clearly, I would ideally have small time steps during the pulses and larger ones inbetween, but how small must these be?! Now the timesteps are 1/100 of the pulse FWHM, which should be sufficient to see something at all, but I do not see anything...

Any further ideas?

[Patrick]

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]

Patrick is right, specifying time points essentially interpolates the solution obtained by the ode solver, it doesn't necessarily modify the steps taken

Answer 2

[Jan]

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

[Joachim]

Hi!

OK. Thank you for your suggestion.

The tolerances may be nonsense. I must admit that I did not think much here.

The function values are not that extreme, however... The pulses have a peak height of ~11, the initial populations are 1 and 0, respectively.

Changing tolerances to ('RelTol', 1e-2,'AbsTol',[1e-4 1e-4]) doesn't make a change. The timespan is 0 to 20000 with features in the driving function as small as 0.1. Scaling the time 10fold did also not change anything.

I boiled the example further down (see below). Now the pulses are several time units in width, height is 1, the function I_drv now enters linearly into the ODE-set (not squared as before). Tolerances are now reasonable.

Nevertheless, with

      options = odeset('RelTol', 1e-2,'AbsTol',[1e-4 1e-4],'maxstep',1);

it works (this time with maxstep 1 instead of 0.1); with

      options = odeset('RelTol', 1e-2,'AbsTol',[1e-4 1e-4]);

it doesn't work, and only the first pulse enters correctly into the ODE solution.

Do I simply have to reduce maxstep and see whether the result changes or not?! Or are there certain rules?!

Perhaps my definition of the pulse-function is bad somehow... But the ploted pulses (using the function) look good.

Best regards, Joachim

function run
%parameters for the driving function
    t_p=10;                     %time of first pulse (ns)
    sigma_p=2;                  %pulse width in ns (sigma) 
    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=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)
            % state N1
        Ndot(1) = -I_drv(t)*0.1*N(1);
            % state N2    
        Ndot(2) = +I_drv(t)*0.1*N(1);

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

%solve ODEs
    options = odeset('RelTol', 1e-2,'AbsTol',[1e-4 1e-4]);%,'maxstep',1);
    [T,Y]=ode15s(@myode,[0 20000],[1 0],options);

    %Plot ODEs    
        f_handle=figure(1);

    subplot(2,1,1)
    title 'driving pulses'
    plot(0:0.001:20000,I_drv(0:0.001:20000),'red');

    subplot(2,1,2)
    title 'state populations'
    plot(T,Y)
    legend('N1', 'N2')
    ylim([0 1])
end

[Joachim]

Any further suggestions, anyone?

[Patrick]

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

[Joachim]

Hi Patrick,

sorry, I must have missed your comment the first time. Yes, I tried calling the solver several times. This works fine, yet is indeed a little unsatisfactory. But it works and is much faster than keeping the small time step for the whole timespan.

Thank you for that comment!

[Richard Brown]

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