vectorizing a for loop in Euler's method (ODE solver)

Phil Isaac — Tue, 11 Aug 2009 02:28:05 -0400

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

Re: vectorizing a for loop in Euler's method (ODE solver)

Bruno Luong — Tue, 11 Aug 2009 07:44:03 -0400

I have my doubt non-linear iterative can be vectorized.

Bruno

MATLAB Central Newsreader - vectorizing a for loop in Euler's method (ODE solver)

vectorizing a for loop in Euler's method (ODE solver)

Re: vectorizing a for loop in Euler's method (ODE solver)