Matlab ode solvers: changing state and specified time

by

Inserting large discontinuities in your ODE in the way you suggest (and the way illustrated by @macduff) can lead to less precision and longer computation times (especially with ode45 – ode15s might be a better option or at least make sure that your absolute and relative tolerances are suitable). You’ve effectively produced a very stiff system. If you want to add some number to the ODEs starting at a specific time keep in mind that the solver only evaluates these equations at specific points in time of its own choosing. (Don’t be mislead by the fact that you can obtain fixed step-size outputs by specifying tspan as more than two elements – all of Matlab’s solvers are variable step-size solvers and choose their true steps based on error criteria.)

A better option is to integrate the system piecewise and append the resultant outputs from each run together:

% t = 0 to t = 10, pass parameter a = 0 to add to ODEs
a = 0;
tspan = [0 10];
[T,Y] = ode45(@(t,y)myfun(t,y,a),tspan,y_0);

% t = 10 to t = 20, pass parameter a = 10 to add to ODEs
a = 10;
[t,y] = ode45(@(t,y)myfun(t,y,a),tspan+T(end),Y(end,:));
T = [T;t(2:end)];
Y = [Y;y(2:end,:)];

% t = 20 to t = 30, pass parameter a = 20 to add to ODEs
a = 20;
[t,y] = ode45(@(t,y)myfun(t,y,a),tspan+T(end),Y(end,:));
T = [T;t(2:end)];
Y = [Y;y(2:end,:)];

The Matlab editor may complain about the array T and Y not being preallocated and/or growing, but it’s fine in this case as they’re growing in large chunks only a few times. Alternatively, if you want fixed output step-sizes, you can do this:

dt = 0.01;
T = 0:dt:30;
Y = zeros(length(T),length(y_0));

% t = 0 to t = 10, pass parameter a = 0 to add to ODEs
a = 0;
[~,Y(1:10/dt+1,:)] = ode45(@(t,y)myfun(t,y,a),T(1:10/dt+1),y_0);

% t = 10 to t = 20, pass parameter a = 10 to add to ODEs
a = 10;
[~,Y(10/dt+1:20/dt+1,:)] = ode45(@(t,y)myfun(t,y,a),T(10/dt+1:20/dt+1),Y(10/dt+1,:));

% t = 20 to t = 30, pass parameter a = 20 to add to ODEs
a = 20;
[~,Y(20/dt+1:end,:)] = ode45(@(t,y)myfun(t,y,a),T(20/dt+1:end),Y(20/dt+1,:));

One could easily convert both of the above blocks of code to more compact for loops if desired.

In both cases your ODE function myfun incorporates the parameter a this way:

function ydot = myfun(t,y,a)
    y(1) = ... % add a however you like
    ...

More Related Contents:

  1. MATLAB: Is it possible to have two event values whilst using ode45?
  2. Absolute error of ODE45 and Runge-Kutta methods compared with analytical solution
  3. Solving an ODE when the function is given as discrete values -matlab-
  4. Matlab ode45. How to change a parameter inside it while calling it?
  5. How to draw a Histogram in Matlab
  6. For loop to split matrix to equal sized sub-matrices
  7. How to normalize a histogram in MATLAB?
  8. Difference between accessing cell elements using curly braces and parentheses
  9. Understanding Matlab FFT example
  10. Nearest-neighbor interpolation algorithm in MATLAB
  11. How to assign values to a MATLAB matrix on the diagonal?
  12. pdist2 equivalent in MATLAB version 7
  13. How do I resize a matrix in MATLAB?
  14. Multi-class classification in libsvm [closed]
  15. Run a script that uses multiple MATLAB sessions
  16. How to create an executable .exe file from a .m file
  17. Analytical Fourier transform vs FFT of functions in Matlab
  18. Index of max and min value in an array
  19. How to deal with overflow and underflow?
  20. What is the closest thing MATLAB has to namespaces?
  21. Cost function in logistic regression gives NaN as a result
  22. Turning y axis upside down in MATLAB
  23. print variable-name in Matlab
  24. Is there a way to see the source of a pcode file (.p) in MATLAB?
  25. Mex files: how to return an already allocated matlab array
  26. MatLab error: cannot open with static TLS
  27. Iterating through struct fieldnames in MATLAB
  28. Save Matlab invisible plot under terminal as an image with same size
  29. Using OpenGL in Matlab to get depth buffer
  30. MATLAB linked list

Leave a Comment