2

3D 関数をベクトル化したいのですが、関数に解析式がありません。たとえば、関数をベクトル化できます

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

のようなことをすることによって

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

そして、スクリプトを実行します

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

ただし、関数に分析式が含まれていないとします。たとえば、

 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

whereODE_fnは、ODE の右辺を吐き出す関数です。したがって、関数は単純に ODE を解くだけなので、関数は明示的にはわかりません。もちろん、ループを使用することもできますが、それらは遅くなります (特にODE を解くために使用する forので、Octave では好まれます)。lsode

のようなものを試す

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

動作しません。また、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
4

1 に答える 1

0

ode45 を使用して微分方程式系を解くことができるため、ODE_fn がベクトル化されている場合は、このアプローチを使用できます。y0 もベクトルである必要があります。

[x1, ..., xn, y1, ..., yn, z1, ..., zn, M1_1-9, ... ,Mn_1-9] である y0 を作成し、x に使用できます。 y、z は適切なインデックス、つまり 1:n、n+1:2*n、2*n+1:3n です。次に、reshape(yin(3*n+1:end),3,3,n) を使用します。しかし、行列の乗算をベクトル化する方法がわかりません。

于 2012-08-09T10:49:53.797 に答える