|
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
|