在Matlab中加速物理模拟计算

Question

我正在做一个用Matlab编写的MR-物理模拟，它模拟布洛赫在一个定义的对象上的方程。对象中的磁化功能每一次更新一次，具有以下功能。

function Mt = evolveMtrans(gamma, delta_B, G, T2, Mt0, delta_t)
    % this function calculates precession and relaxation of the
    % transversal component, Mt, of M
    delta_phi = gamma*(delta_B + G)*delta_t;        
    Mt = Mt0 .* exp(-delta_t*1./T2 - 1i*delta_phi);
end  

这个函数在整个代码中只占很小的一部分，但是被调用了多达250.000次，从而减慢了代码和整个模拟的性能。我曾想过如何加快计算，但还没有想出一个好的解决办法。有一条线路非常耗时，约占整个模拟时间的50% - 60%。这是台词，

Mt = Mt0 .* exp(-delta_t*1./T2 - 1i*delta_phi);

哪里

Mt0 = 512x512矩阵

delta_t =标量

T2 = 512x512矩阵

delta_phi = 512x512矩阵

如能提出加快计算速度的建议，我将不胜感激。

下面的更多信息，

在仿真过程中，每个时间步骤都调用函数evovleMtrans。用于调用函数的参数是，

gamma =常数。(旋转磁常数)

delta_B =磁场值

G =梯度强度

T2 =具有T2值的512x512矩阵

Mstart.r =一个值为M.r的512x512矩阵有最后一个时间步骤

delta_t =自上一次计算M.r以来具有时间差的标量

在模拟过程中，这些参数的唯一变化是：G、Mstart.r和delta_t。其余的在模拟过程中不会改变它们的值。

下面的部分是调用函数的主代码中的部分。

        % update phase and relaxation to calcTime

        delta_t = calcTime - Mstart_t;
        delta_B = (d-d0)*B0;
        G = Sq.Gx*Sq.xGxref + Sq.Gz*Sq.zGzref;

        % Precession around B0 (z-axis) and B1 (+-x-axis or +-y-axis)
        % is defined clock-wise in a right hand system x, y, z and
        % x', y', z (see the Bloch equation, Bloch 1946 and Levitt
        % 1997). The x-axis has angle zero and the y-axis has angle 90.
        % For flipping/precession around B1 in the xy-plane, z-axis has
        % angle zero.
        % For testing of precession direction:
        % delta_phi = gamma*((ones(size(d)))*1e-6*B0)*delta_t;
        M.r = evolveMtrans(gamma, delta_B, G, T2, Mstart.r, delta_t);
        M.l = evolveMlong(T1, M0.l, Mstart.l, delta_t); 

Answer 1

这并不令人惊讶。

“单线”是一个矩阵方程。它实际上是1,024个联立方程。

按照Jannick的说法，第一个词意味着按元素进行除法，所以“δ_t/Ti，j”。矩阵乘以标量是O(N^2)。矩阵加法为O(N^2)。矩阵的指数将是O(N^2)。

我不确定我是否也在里面看到了一个复杂的论点。这是否意味着具有真实和想象条目的复杂矩阵？你的方程是否简化为实部和虚部？这意味着计算次数的两倍。

你最好的希望是尽可能利用对称性。如果你所有的矩阵都是对称的，你削减你的计算差不多一半。

如果可以的话，使用并行化。

算法的选择也会产生很大的影响。如果您使用的是显式欧拉集成，则由于稳定性问题，您可能有时间步长限制。这就是你有25万步的原因吗？也许通过一个更稳定的集成模式，一个更大的时间步骤是可能的。考虑一种具有误差修正的高阶自适应方案，如5阶Runge Kutta。

Answer 2

有几种可能提高代码的速度，但我所看到的都有一个警告。

数值ode积分

第一种可能是通过数值微分方程解来改变你的解析解。这有几个优点

1. 解析解包含复指数函数，计算费用高，而微分方程只包含乘法和加法。(d/dt u= -a u => u=exp(-at))
2. 有许多内置的解决方法，为matlab可用，他们通常是相当快(例如ode45)。然而，内置的所有步骤都使用可变的步长。这提高了速度和准确性，但如果你真的需要一个固定的、等距的时间点网格，那将是一个问题。这里是非官方的固定步长求解器。

首先，您还可以尝试只使用euler步骤，方法是替换

M.r = evolveMtrans(gamma, delta_B, G, T2, Mstart.r, delta_t);

通过

delta_phi = gamma*(delta_B + G)*t_step;
M.r += M.r .* (1-t_step*1./T2 - 1i*delta_phi);

然后，您可以通过预先计算所有常量值(例如one_over_T1=1/T1 )，将delta_phi移出循环，从而进一步改进这一点。

警告:你必须有一个最小的步长，否则精度就会受到影响。因此，这只是一个好主意，如果你的时间间隔是相当好的。

减少时间点

你应该认真分析你是否真的需要这么多时间点。在我看来，你需要这么多要点，这似乎有点令人费解。正如您所知道的完整的分析解决方案，您可以自由选择如何采样时间，并可能利用这一点，以您的优势。

去fortran

这看起来似乎是一大步，但在我的经验中基本(简单的循环，矩阵运算等)matlab代码可以相对容易地翻译成fortran逐行.除了我的第一点外，这将特别有帮助。如果您仍然想使用完整的解析解，可能这里没有什么收获，因为exp在matlab中已经非常快了。