-1

I have rather lengthy system of equations saved in a function file, system3. I need one of the parameters in this system to change dependant on time. I have created a second function, calculate_a1, that produces a vector of the parameter a1 at each of my 401 time points. 

  tResult = []; 
  xResult = []; 
  tStep = linspace(0,400,401); 
  y0    = [IC]; 
  alpha = calculate_a1();

  for index = 2:numel(tStep)

  % Integrate:

  a1 = alpha(1,index);
  t    = tStep(index-1:index);
  sol = ode45(@system3,t,y0,a1)

  % Collect the results:

  tResult = cat(1, tResult, t);
  xResult = cat(1, xResult, x);

  % Final value of x is initial value for next step:

  y0 = x(end, :);

  end

Up until the sol line, this works fine, but I'm struggling to export a1 with ode45 so that it can be used to solve system3. Any help would be appreciated.

horchler
  • 18,384
  • 4
  • 37
  • 73
Avg
  • 21
  • 5
  • Questions seeking debugging help ("why isn't this code working?") must include the desired behavior, a specific problem or error and **the shortest code necessary to reproduce it in the question itself**. Questions without a clear problem statement are not useful to other readers. – sco1 Jan 05 '16 at 15:37
  • you know that the RHS (in your case system3) can accept time as a parameter right? Why can't the time-dependence be implemented inside the system3 function? http://se.mathworks.com/help/matlab/ref/ode45.html?refresh=true – alexandre iolov Jan 05 '16 at 16:26
  • "I need one of the parameters in this system to change dependant on time." – if it depends on the independent variable, time, then by definition it's not a parameter. As mentioned by @LutzL, you may be able to use interpolation inside your integration function – see [this question/answer for an example](http://stackoverflow.com/q/19732442/2278029). If the "parameter" only changes at a few key points in time then you can try integrating each part piecewise – see [here](http://stackoverflow.com/q/17356037/2278029) or [here](http://stackoverflow.com/a/17370733/2278029). – horchler Jan 05 '16 at 22:39

1 Answers1

1

That will not work. ode45 will need to evaluate this parameter at a lot more time points. First, the time list you give only determines the output, its values are interpolated from the internal samples which are situated at dynamically regulated sample time points. And second, every internal time step consists of 5 evaluations of the ODE function at different fractions of the time step.

The best you can do is to give the function behind calculate_a1 as parameter, so that in every evaluation of the ODE function the correct a1 corresponding to this exact time can be calculated.

Interpolation is not that good an idea since at the nodes the step size regulator will "feel" the limited smoothness as stiffness of the system and thus regulate the step size down, which greatly increases computation times and the accumulation of floating point errors.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • How would I go about giving calculate_a1 as a parameter? The value of a1 isn't calculated, it comes from data that I have collected. At the moment, I have a vector the same length as tStep that has all of the values for a1 in it. Apologies if i'm going about asking these questions in the wrong way, I'm new to the worlds of both MatLab and StackExchange! – Avg Jan 06 '16 at 19:38
  • Then you need to go via interpolation and accept the inefficiencies or you hand-code a multi-step method for the discretization of `tStep`, I do not remember seeing library functions for that. – Lutz Lehmann Jan 06 '16 at 19:53
  • Running the code using interpolation, I get the error `Unable to meet integration tolerances without reducing the step size below the smallest value allowed` and an empty plot! . Is this likely to be from the use of interpolation or an error in my system of ODE's? – Avg Jan 06 '16 at 20:37