CODEWITHSUNDEEP
matlaboctavevectorization
db1234
db1234

Reputation: 807

vectorizing a non-explicit function in Matlab and Octave

I want to vectorize a 3D function, but the function does not have an analytical expression. For instance, I can vectorize the function

 F(x, y, z) = (sin(y)*x, z*y, x*y)

by doing something like 

 function out = Vect_fn(x, y,z)
     out(1) = x.*sin(y);
     out(2) = z.*y;
     out(3) = x.*y;
 end

And then running the script

 a = linspace(0,1,10);
 [xx, yy, zz] = meshgrid(a, a, a);
 D = Vect_fn(xx, yy, zz)

However, suppose the function does not have an analytical expression, for example

 function y = Vect_Nexplicit(y0)
      %%%%%%y0 is a 3x1 vector%%%%%%%%%%%%%%
      t0 = 0.0; 
      tf = 3.0;
      [t, z] = ode45('ODE_fn', [t0,tf], y0);
      sz = size(z);
      n = sz(1);
      y = z(n, :);
 end

where ODE_fn is just some function that spits out the right-hand side of a an ODE. Thus the function simply solves an ODE and so the function is not known explicitly. Of course I can use a for loop, but those are slower (esp. in Octave, which I prefer since it has lsode for solving ODEs)

Trying something like

a = linspace(0,1,10);
 [xx, yy, zz] = meshgrid(a, a, a);
 D = Vect_Nexplicit(xx, yy, zz)

does not work. Also here is the code for ODF_fn:

 function ydot = ODE_fn(t, yin)

 A = sqrt(3.0);
 B = sqrt(2.0);
 C = 1.0;

 x = yin(1, 1);
 y = yin(2,1);
 z = yin(3, 1);

 M = reshape(yin(4:12), 3, 3);

 ydot(1,1) = A*sin(yin(3)) + C*cos(yin(2));
 ydot(2,1) = B*sin(yin(1)) + A*cos(yin(3));
 ydot(3,1) = C*sin(yin(2)) + B*cos(yin(1));

 DV = [0 -C*sin(y) A*cos(z); B*cos(x) 0 -A*sin(z); -B*sin(x)      C*cos(y) 0];

 Mdot = DV*M;

 ydot(4:12,1) = reshape(Mdot, 9, 1);


 end

Upvotes: 2

Views: 270

Answers (1)

simmmons
simmmons

Reputation: 211

You can solve systems of differential equations with ode45, so if ODE_fn is vectorised you can use your approach. y0 just needs to be a vector too.

You can create a y0 that is [x1, ..., xn, y1, ..., yn, z1, ..., zn, M1_1-9, ... ,Mn_1-9] and then use for x, y, z just the appropriate indexes i.e. 1:n, n+1:2*n, 2*n+1:3n. Then use reshape(yin(3*n+1:end),3,3,n). But I am not sure how to vectorize the matrix multiplication.

Upvotes: 0

Related Questions