Help with modifying built-in ode45

Hi I need to solve some ODEs using ode45 function. But, I want to insert a little modification. I get the original code by typing edit ode45. Then I inserted my modification and saved it as a different function and called it to solve my odes. But it didn't work. I am new to MATLAB. Can anyone help me to properly insert the modification? They are just couple of if statements limiting couple of outputs bellow zero value. Bellow is a custom RK-4 code with my modification that works, but it's not adaptive like ode45. I want to apply the same exact modification to original ode45 function. Please somebody help.

function [tout, yout] = cnoderklimint(FCLsim, tini, tfin, inic)
h = 0.0001; % h is timestep
time = tini;
tout(1,1) = time;
y = inic;
yout(1,:) = inic;
yp = feval(FCLsim,time,inic); % k1
step = 0;
while(time<tfin)
step = step + 1;

time = time + 0.5 * h;

ytem = y + 0.5 * h * yp;

ak1= feval(FCLsim,time,ytem);    % k2

ytem = y + 0.5 * h * ak1;

ak2 = feval(FCLsim,time,ytem);   % k3

time = time + 0.5 * h;

ytem = y + h * ak2;

ak3 = feval(FCLsim,time,ytem);   % k4

yout(step+1,:) = y+h*(yp+2*(ak1+ak2)+ak3) / 6.0;

y = (yout(step+1,:))';
%%%% Modification starts
if y(3,1)<0

      y(3,1)=0;

end

if y(4,1)<0

      y(4,1)=0;

end

if y(7,1)<0

      y(7,1)=0;

end

if y(8,1)<0

      y(8,1)=0;

end
%Modification ends
yp = feval(FCLsim,time,y);

tout(step+1,1) = time;
end

[Torsten]

Don't modify built-in MATLAB functions, espcially if you are new to MATLAB.

If you want to limit certain variables, do this in your function where you define the ODEs, not in the ODE45 source code.

Best wishes

Torsten.

[Jan]

Please do not post a question multiple times. See https://www.mathworks.com/matlabcentral/answers/335872-help-with-modifying-built-in-ode45. Delete one of the threads before somebody adds an answer.

@Torsten: The OP saved the modified function with a different file name. Then the original function is not destroyed. This is fine even for beginners.

[Nick]

If I do need to modify a function, what is the best practice for this? For example, if I want to modify ode45 then I would save a copy with a different name and modify that. But I get errors because the internal functions like odearguments cannot be found. Should I be saving the modified function in the same folder where ode45 is kept?

[Steven Lord]

IMO the best practice is not to modify MathWorks provided functions. Doing so could impact any code (yours or other MathWorks provided functions) that call the function you modified. Those types of changes in behavior can be very difficult to diagnose (especially if you modified the function months or years ago and forgot about your changes.) And if you've actually modified the files included in the product, installing an Update release could overwrite your modifications if it includes a bug fix that impacts that file.

Consider instead pre-processing the inputs to the MathWorks provided functions, post-processing the outputs from those functions, or (in some cases) asking here on MATLAB Answers if there's functionality in the MathWorks provided functions that you're unaware of that could help you achieve your goal.

For example, in the original question there's an option the poster could use to keep the components of the solution non-negative (which I described in an answer almost 7 years ago.)

[Nick]

Ok. What I am wanting to do is terminate integration after a certain number of failed steps or function evaluations

[Torsten]

Ok. What I am wanting to do is terminate integration after a certain number of failed steps or function evaluations

Be sure: The solver will do that for you early enough.

[Nick]

Oh is there functionality for that?

[Torsten]

If the stepsize in the independent variable becomes too small because the absolute and/or relative error you prescribe in AbsTol and RelTol cannot be met or NaN or Inf values are returned from you ODE function, the solver will automatically quit.

[Nick]

The problem I am facing is a chattering in my ODEs where that mechanism doesn't stop integration. It seems that the tolerance is always able to be satisfied over this chattering and so the step sizes are consistently really small and progress becomes extremely slow

[Steven Lord]

Have you tried a stiff ODE solver?

[Walter Roberson]

It sounds as if you have a stiff system but you are not using a stiff solver.

[Nick]

The solver I usually use is ode113 because it is the fastest for my problem except for when these occasional cases where the integration slows down almost to a stop. I just tried one of the those cases with ode15s and it returns the warning that the tolerances cannot be met. Does that mean it is stiff? Or no since it wasn't able to solve it?

[Walter Roberson]

Experiment with ode23s

[Torsten]

I just tried one of the those cases with ode15s and it returns the warning that the tolerances cannot be met. Does that mean it is stiff?

You use a stiff solver - thus stiffness cannot be the reason. Your equations seem to be inherently difficult to solve in some cases. But now it seems you found the correct solver that quits when things become hopeless.

[Paul]

Potentially relevant blog post just showed up today.

Answer 1

[Jan]

Limiting the trajectory to certain values does not happen inside the integrator. You do not have to modify ode45.m, but the function to be integrated, in your case FCLsim. But then this function is not smooth anymore. As many otehr integrators, Matlabs ODE-suite cannot handle discontinuities and the step size controller might either stop the intergation or even worse proceed and reply an inaccurate result. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. Therefore: Don't do this. Call the integrator in a loop and use an event function to modify a parameter provided as argument of an anonymous function (see http://www.mathworks.com/matlabcentral/answers/1971) .

Answer 2

[Steven Lord]

I strongly discourage you from modifying functions included in MATLAB -- and in this case, it's not necessary for you to do so.

Call odeset to create an options structure with the 'NonNegative' option set to the indices of the components that you want to keep non-negative (in this case [3,4,7,8]) and pass that options structure into ode45 as the fourth input argument.