Issue with fractional step-sizes

4 次查看(过去 30 天)
Patrick
Patrick 2011-3-1
I am writing a m.file with a series of numerical methods which can be called. Several of the methods (Runge-Kutta-Merson and Richardson Extrapolation) have fractional step-sizes. When I run the program I receive the following error "??? Subscript indices must either be real positive integers or logicals". I have tried sizing the step sizes up to avoid the fractions but in some cases it is not possible and I get errors. Part of the for-loop looks like this:
y1=y0*ones(n,1);
y2=y0*ones(n,1);
y=y0*ones(n,1);
for i=1:n-1
y1(i+1)=y(i)+h*dydt(t(i),y(i),varargin{:});
y2(i+1/2)=y(i)+(h/2)*dydt(t(i),y(i),varargin{:});
y2(i+1)=y2(i+1/2)+(h/2)*dydt(t(i+1/2),y(i+1/2),varargin{:});
y(i+1)=2*y2(i+1)-y1(i+1);
end
Any help is greatly appreciated!

回答(6 个)

Matt Tearle
Matt Tearle 2011-3-1
You're mixing up the "physical" steps and the array indexing. In RK methods, the intermediate steps are evaluated between collocation points, but that doesn't mean you need to use fractional indices to the arrays. y(k) stores the kth step of y, which happens to be located at t = t(k) = t0 + k*h; (or maybe (k-1), depending on whether you're storing t0 as t(k=1)).
You can evaluate the rate equations at t = t(k) + h/2, y = y(k) + (h/2)*(some slope). But again, that's not the same thing as t = t(k+1/2), y = y(k+...).
As an aside: do you really need to keep all the intermediate calculations at each step? Typically these are stored as temporary variables that get overwritten at each step. That avoids a lot of indexing issues.

Sean de Wolski
Sean de Wolski 2011-3-1
How can you store data in the i.5th location?

Brett Shoelson
Brett Shoelson 2011-3-1
y(1) refers to the first element of y. y(2) is its second element. What does y(1.5) refer to?
Your problem is in your use of y(i+1/2) and t(i+1/2). You need to rework your code so that your variable indices are integers. For instance:
incr = 0;
for ii = 1:0.5:n
incr = incr + 1;
y1(incr) = f(ii);
end
Cheers,
Brett

Oleg Komarov
Oleg Komarov 2011-3-1
You have to keep separate the concept of your step size from the position of an element in an array:
A = [.25 .5 .75 1]
pos 1 2 3 4
If n say 100, then my A would be 100*4, therefore I'll have 400 positions.
Oleg

Patrick
Patrick 2011-3-1
I understand what you are getting at, and its where I thought I was going wrong. Matt, I have been trying to do what you suggested, but how would I choose the slope to multiply the (h/2) by?

Patrick
Patrick 2011-3-1
Actually, I think i just got it.
  1 个评论
Matt Tearle
Matt Tearle 2011-3-1
OK. But either way, the answer is "from the specific RK algorithm". If you don't need to keep all the intermediate steps, you can often just use temp variables.
slope1 = odefun(t(k),y(k))
ytemp = y(k) + (h/3)*slope1;
slope2 = odefun(t(k)+h/3,ytemp);
ytemp = y(k) + (h/2)*slope2;
etc
y(k+1) = y(k) + (h/42)*(17*slope1 - pi*slope2 + etc);

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Loops and Conditional Statements 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by