【问题标题】:vectorizing a non-explicit function in Matlab and Octave在 Matlab 和 Octave 中对非显式函数进行矢量化
【发布时间】:2012-07-27 13:37:31
【问题描述】:

我想对一个 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

其中ODE_fn 只是一些输出 ODE 右侧的函数。因此,该函数仅求解 ODE,因此该函数不明确。当然,我可以使用for 循环,但速度较慢(尤其是在 Octave 中,我更喜欢它,因为它有 lsode 用于求解 ODE)

尝试类似的东西

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

【问题讨论】:

  • 如果没有解析解,并且您必须为每个新输入求解一个微分方程,我看不出通过某种循环构造如何实现。 (matlab中的arrayfun等本质上也是一个循环)
  • 我认为如果 ODE_fn 被正确矢量化,您当前的方法可能会奏效。

标签: matlab octave vectorization


【解决方案1】:

您可以使用 ode45 求解微分方程组,因此如果 ODE_fn 被矢量化,您可以使用您的方法。 y0 也需要是一个向量。

您可以创建一个 y0,即 [x1, ..., xn, y1, ..., yn, z1, ..., zn, M1_1-9, ... ,Mn_1-9] 然后使用对于 x, y, z 只是适当的索引,即 1:n、n+1:2*n、2*n+1:3n。然后使用 reshape(yin(3*n+1:end),3,3,n)。但我不确定如何向量化矩阵乘法。

【讨论】:

  • 谢谢。将我的特定形式的 ODE_fn 矢量化可能有点困难,但我可以尝试。也很高兴了解 t0 和 tf。
  • 只要把它看作一个解耦微分方程系统。
  • 解耦不是正确的词。对于二维方程,非耦合均值 (xdot, ydot) = (f(x), g(y))。请注意,系统 (xdot, ydot) = (f(x, y),g(x, y)) 不是解耦的。因为我正在进化一个矩阵,所以我所拥有的更加复杂......
  • 我在原帖中添加了 ODE_fn
  • 我刚刚意识到不是 t0 和 tf 需要是向量,而是 y0。你应该创建一个解耦的方程组。我会相应地改变我的答案。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2011-07-08
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多