for loop ode 45 changing parameter values

I'm trying to run this code and use only the last value of the tspan for the ode45, however I always get 2 values for the final value of the tspan (when t=100)

This is the code:

ABP=linspace(40,170,131);
ti=0; tf=100;
 %change in ABP at every loop   

       for i=1:1:length(ABP)    
           sol = ode45(@first,ti:1:tf, [0 0 0 0 0.205 62.516 12 13.201 0.2 1],[],ABP(i)); %(function, timespan, initial condition for xq,xc,xm1,xm,Ca,P1,V_sa,P2)
       end

%%========= FUNCTION ==================
function dvdt = first(t,v,ABP)
%code 
...ABP used in a set of differential eq...
   if t==100
van= (Vsa+Vla)/(Vsab+Vla) 
fileID=fopen('van.txt','a');
fprintf(fileID,' %4.3f\n',van);
fclose(fileID);

   end
end

Hence, I should end up having a .txt file with 131 values of van, one for each value of ABP. Instead, I end up with double the amount. I also printed the ABP values in the txt file to check, and indeed there were two different (very similar) values of van for each value of ABP. Any suggestions on how to fix this?

[Jan]

I still cannot guess, what "end up with double the amount" means. You get "2 values" of what?

[gorilla3]

Instead of getting 131 values for "van" I get 262. i.e. I get 2 results for each value of ABP

Answer 1

[David Goodmanson]

Hi gorilla3,

Here is the simplest possible example of what ode45 is doing, that also inputs an extra parameter. To find the times that ode45 employs, I used the low-rent method of uncommenting the commands in ‘first’ and letting the results run out in the command window. If you do that you will see that the last time occurs twice, although it of course only occurs once in the output sol.x.

I suppose that Vsa, Vsab, etc. depend on v, APB, and possibly t. I think you will probably have to take some code outside of ode45 (a copy if some of it was actually needed to help get dvdt) and do the calculation within the ABP loop, using sol.y(end) and (if time is needed) sol.x(end).

ABP = .1;
sol = ode45(@first,[0:100],[1],[],ABP);
plot(sol.x,sol.y,'o-')

function dvdt = first(t,v,ABP)
t
dvdt = -v*ABP
end

[gorilla3]

Thanks for your reply. However, I don't understand the solution you suggest... in your code you removed the for loop that I need to change values of ABP every time. Also, I can't perform some calculations outside the function because all other parameters are internally stored in the function itself.

Maybe a solution would be finding another way to change ABP values at every loop. Any suggestions?

[David Goodmanson]

I wanted to show was going wrong with ode45 for demonstration purposes, and there was no need to keep the for loop and repeat that 131 times.

Here is what I am suggesting in more detail. Within your function 'first' you are doing some calculation to find van. That calculation depends on v and (maybe) t. If you had v and t, you could calculate van just as well outside of 'first' as inside.

You want van for the last values of t and v that were calculated by ode45. That information is saved in sol, so you can do the calculation outside of 'first' but still inside the loop.

Also, you don’t have to define the ‘first’ function every time in the loop. You can save it as a function mfile somewhere on the Matlab path where Matlab can find it, or if you are doing this in a script you can put the function at the bottom, making sure it has and ‘end’ statement at the end.

   ABP=linspace(40,170,131);
   ti=0; tf=100;

 for i=1:1:length(ABP)    %change in ABP at every loop  

   % (function, timespan, initial condition for xq,xc,xm1,xm,Ca,P1,V_sa,P2)
  sol = ode45(@first,ti:1:tf,[0 0 0 0 0.205 62.516 12 13.201 0.2 1],[],ABP(i));
   vfinal = sol.y(end);
   tfinal = sol.x(end);  
   van = < calculation involving vfinal, tfinal, constants >  

   fileID=fopen('van.txt','a');
   fprintf(fileID,' %4.3f\n',van);
   fclose(fileID);
end

========= FUNCTION ==================
function dvdt = first(t,v,ABP)
%code 
end

[gorilla3]

Thank you very much for elaborating your reply. It's very helpful!

Answer 2

[Jan]

Fo not include an output inside the function to be integrated. Remember that the integrator can call this function with different inputs to evaluate one step. The integrator can even reject steps.

Use the OutputFcn instead, which is called after a successful step:

options = odeset('OutputFcn', @yourOutput)

This has been suggested in your other thread https://www.mathworks.com/matlabcentral/answers/371408-selecting-last-value-of-tspan-ode45#answer_294967 already.

For completeness I repeat, that providing the parameters as input of ODE45 is deprecated for over 17 years now, see http://www.mathworks.com/matlabcentral/answers/1971> .