MATLAB Answers

Problem with ode15s / ode45

79 views (last 30 days)
Joachim
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

  0 Comments

Sign in to comment.

Accepted Answer

Star Strider
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.

  3 Comments

Joachim
Joachim on 30 Jul 2012
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
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
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

Sign in to comment.

More Answers (1)

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

  5 Comments

Show 2 older comments
Patrick
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?
Joachim
Joachim on 10 Aug 2012
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
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

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!