for loop ode 45 changing parameter values

4 次查看（过去 30 天）
gorilla3
gorilla3 2017-12-13
评论： gorilla3 2017-12-15
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?
  2 个评论
Jan
Jan 2017-12-14
I still cannot guess, what "end up with double the amount" means. You get "2 values" of what?
gorilla3
gorilla3 2017-12-15
Instead of getting 131 values for "van" I get 262. i.e. I get 2 results for each value of ABP

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

David Goodmanson
David Goodmanson 2017-12-13
编辑：David Goodmanson 2017-12-13
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
  3 个评论
显示 1更早的评论
David Goodmanson
David Goodmanson 2017-12-15
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
gorilla3 2017-12-15
Thank you very much for elaborating your reply. It's very helpful!

请先登录，再进行评论。

更多回答（1 个）

Jan
Jan 2017-12-14
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> .

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by