loop for with ode45

Question

tony 2013-1-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/59695-loop-for-with-ode45

在 MATLAB Online 中打开

hi , i've got a problem..

I want to solve a EDO , so i used ODE45, i've got two function :

*function dy=f_iso(t,y)

global Q b K n E sigma0 epoint

if abs(y(1))-Q*(1-exp(-b*y(2)))-sigma0<0

dy(2)=0

else

dy(2)=((abs(y(1))-Q*(1-exp(-b*y(2)))-sigma0)/K).^n

end

dy(3)=sign(y(1))*dy(2);

dy(1)=E*(epoint-dy(3));

dy=dy';

end*

And my main function :

[Tn,N]=ode45('f_iso',[0 10],[0 0 0]);

sigma1=N(:,1);

ep1=N(:,3);

e1=ep1+sigma1./E;

It's the first court , then i did the second : [Th,h]=ode45('f_iso',[10 30],[N(89,:]);

sigma2=h(:,1);

ep2=h(:,3);

e2=ep2+sigma2./E;

plot(e1,sigma1,e2,sigma2)

and i want to do it 10 times , so i did a loop for :

k=3

for i=1:k

 Ti=10
    Tf=30
    epoint=(-1)^i*0.001;
    C=N(length(N),:)
    [T,b]=ode45('f_iso',[Ti Tf],[C])
    Ti=Ti+20
    Tf=Tf+20
    C=b(length(b),:);

end

but i doest work and i don't know how to use plot with a loop for...

tranks :)

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2013-1-24

Please format your code properly. Currently it is unnecessarily hard to read your question.

You did not answer my questions for clarifications in your related question http://www.mathworks.com/matlabcentral/answers/59582-ode45-errer-input-argument-y-is-undefined. There I've explained already, that ODE45 is not sufficient to solve a discontinous system. I strongly suggest to consider this.

Please explain what "it does not work" exactly means: Do you get an error message or does the result differ from your expectations.

请先登录，再进行评论。

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

Answer 1

Jan 2013-1-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/59695-loop-for-with-ode45#answer_72250

编辑：Jan 2013-1-24

在 MATLAB Online 中打开

You reset Ti, Tf and C in each iteration to the same values. Perhaps you want this:

global epoint   % !!!
Ti = 10;
Tf = 30;
C  = N(size(N, 1), :);
for i=1:k
  epoint = (-1)^i*0.001;   % Not used at all?
  [T, b] = ode45(@f_iso, [Ti Tf], C);
  Ti = Tf
  Tf = Ti+20
  C  = b(size(b, 1), :);
end