Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: vectorizing a for loop in Euler's method (ODE solver)
Date: Tue, 11 Aug 2009 02:28:05 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 33
Message-ID: <h5qkvl$92r$1@fred.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1249957685 9307 172.30.248.35 (11 Aug 2009 02:28:05 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 11 Aug 2009 02:28:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1930752
Xref: news.mathworks.com comp.soft-sys.matlab:562280


As an educational exercise, I'd like to know if the for loop in the following M-file code can be vectorized:

%%%
t0=0;
y0=4;
tMax = 5;
steps= 20;
rhsFunc = @(t,y)(t.* y.^2 .* sin(2*y.*t));
y=zeros(1,steps+1);
y(1) = y0;
deltaT = (tMax-t0)/steps;
t = t0:deltaT:tMax; 
for i = 1:steps
  y(i+1) = y(i) + rhsFunc(t(i),y(i))*deltaT;
end
plot(t,y);
%%%

This is an implementation of Euler's method, a simple first order ODE solver.

I understand how to vectorize a loop using "filter" (or in some cases "cumsum" or "cumprod") that creates a vector whose elements depend on the previous element, in the case where those elements depend *linearly* on the previous elements. That is explained clearly in Support - Code Vectorization Guide:

http://www.mathworks.com/support/tech-notes/1100/1109.shtml

In the example code that I have provided, I have purposely made "rhsFunc" a non-linear function of y. I am having trouble working out how to vectorize that loop in this case, or indeed in the general case where "rhsFunc" is any (smooth) non-linear function of two variables, in particular the variable y.

Some additional comments:
- I know about the ODE solvers ode23, ode45, etc. I'm more interested in how to vectorize, rather than the context.
- The code I have provided already works fine. I have used a small number of steps, but if you change to "steps = 20000;", the plot produced is not bad (compare with ode45) and still runs in good time (of course not as efficient as Runge-Kutta!), and so there may be no point in vectorizing at all. Again, I'm interested in the question of how to vectorize, rather than the context.

Cheers,

Phil