Numerical Simulation of Bloch Equations for
Dynamic Magnetic Resonance Imaging
Arijit Hazra∗, Gert Lube∗, Hans-Georg Raumer

arXiv:1703.05985v3 [math.NA] 10 Sep 2017

Institute for Numerical and Applied Mathematics,
Georg-August-University Göttingen, D-37083, Germany

Abstract
Magnetic Resonance Imaging (MRI) is a widely applied non-invasive imaging modality
based on non-ionizing radiation which gives excellent images and soft tissue contrast of
living tissues. We consider the modified Bloch problem as a model of MRI for flowing spins in an incompressible flow field. After establishing the well-posedness of the
corresponding evolution problem, we analyze its spatial semi-discretization using discontinuous Galerkin methods. The high frequency time evolution requires a proper explicit
and adaptive temporal discretization. The applicability of the approach is shown for
basic examples.
Keywords: Magnetic resonance imaging, Bloch model, FLASH-technology, flowing
spins, incompressible medium, discontinuous Galerkin method, explicit Runge-Kutta
methods

1. Bloch Model for Magnetic Resonance Imaging
Magnetic Resonance Imaging (MRI) is a non-invasive imaging modality based on nonionizing radiation [1]. It gives excellent images and soft tissue contrast of living tissues.
During the experiment, the object to be studied is placed in a static magnetic field of
high strength B0 . This induces a macroscopic nuclear magnetization M0 in the direction
of magnetic field, known as the equilibrium magnetization. The direction of equilibrium
magnetization is called the longitudinal direction, generally denoted by the z-axis as in
Fig. 1.
In order to get a response from the object, the equilibrium magnetization is perturbed
by applying a short radio-frequency (RF) pulse B1 in the transverse plane (denoted by
the xy-plane in Fig. 1) with the excitation carrier frequency of the RF pulse equal to
the Larmor frequency of protons. As a result, the equilibrium magnetization is flipped
from its initial position towards the transverse plane. The perturbation of equilibrium

∗ Corresponding

authors
Email addresses: arijit.hazra@stud.uni-goettingen.de (Arijit Hazra),
lube@math.uni-goettingen.de (Gert Lube), hansgeorg.raumer@stud.uni-goettingen.de
(Hans-Georg Raumer)
Preprint submitted to Applied Numerical Mathematics
September 12, 2017

Equilibrium

Excitation

Precession

Decay

+

Recovery

Figure 1: Schematic of a pulsed MR experiment on net Magnetization M. Left: M align along the static
magnetic field B0 . Second from left: RF excitation perturbs the magnetization from the longitudinal
direction. Middle: Precession of magnetization as a result of, Second from right: T2 decay. Right: T1
recovery.

magnetization depends on the duration and magnitude of the RF field. After the RF
pulse is switched off, the magnetization precesses towards its equilibrium position. This
process is known as relaxation. The relaxation of magnetization is governed by two time
constants - spin-lattice time relaxation T1 , spin-spin time relaxation T2 , as illustrated in
Fig. 1.
On top of that magnetic field gradients need to be applied to obtain cross-sectional
images. The location and thickness of the slice is determined by applying a slice-selection
gradient Gz . After that, the object is spatially encoded via the application of additional
magnetic field gradients Gx , Gy in the transversal directions. The magnitude and time
of application of magnetic gradients G depend on the experimental requirements, see
Sec. 5 for the description of a typical pulse sequence.
The Bloch equation [2] for MRI combines all of the above components: static magnetic field B0 , time-dependent magnetic gradients G(t) and the RF pulse B1 (t) =
(Bx (t), By (t), 0)T and is given by
Mx
My
M0 − Mz
dM
êz −
êx −
êy ,
= γM × B +
dt
T1
T2
T2
B(t, r) = (Bx (t), By (t), Bz (t, r))T , Bz (t, r) := B0 + G(t) · r.

(1.1)

Emitted energy due to the precession of the magnetizations is converted into an electric
signal in the receiver coil of an MRI system which is manipulated further for image
reconstruction. The received signal can be expressed as
Z
Rt
Mxy (r, t)e−iγB0 t e−iγ 0 (xGx +yGy )dτ dr
sr (t) ∝
(1.2)
Ω

where Mxy = Mx + iMy . The received signals are typically demodulated in frequency
by γB0 using phase-sensitive detection before being used for image reconstruction. The
resultant signal expression after demodulation is sd (t) = sr (t)eiγB0 t [1].
The demodulated signal corresponds to the solution of the Bloch equation in a frame
rotating clockwise about z-axis with an angular frequency γB0 , i.e, the Larmor frequency
of the isocenter of the object. The Bloch equation in the rotating frame conceptually
simplifies the RF excitation effect in MRI and eliminates the static B0 field from the

2

expression of the external magnetic field and is given by
dM′
Mx′ êx′ + My′ êy′
(M0 − Mz )êz
−
= γM′ × Beff +
dt
T1
T2


  
1
γBz′
−γBy′
−
0
Mx′
 T2

 
− T12
γBx′  My′  +  0  ,
= −γBz′
M0
Mz
− T11
γBy′ −γBx′
T1

(1.3a)

(1.3b)

where M′ (t, r) = (Mx′ , My′ , Mz )T (t, r) represents the magnetization in the rotating
frame and the effective magnetic field is given by Beff (t, r) = (Bx′ , By′ , Bz′ )T (t, r);
Bz′ (t, r) := G(t) · r.
In order to study the effect of fluid flow in an MRI experiment, the transport of
magnetizations due to flow field u(t, r) must be taken into account and is modeled by
the modified Bloch equation [3]:
Mx êx + My êy
(M0 − Mz )êz
∂M
−
.
+ (u · ∇)M = γM × B +
∂t
T1
T2

(1.4)

There will be no diffusion term as effect of diffusion is negligible, but see Sec. 2.
Due to signal demodulation, the signal acquired from the flowing object in an MRI
experiment is equivalent to an equation where magnetizations and the magnetic field
(Beff ) are written in the rotating frame with velocity u(t, r) kept in laboratory frame [4]:
Mx′ êx′ + My′ êy′
(M0 − Mz )êz
∂ M′
−
.
+ (u · ∇)M′ = γM′ × Beff +
∂t
T1
T2

(1.5)

For notational simplicity the ′ will be omitted in later sections. Hereafter, magnetization
M will always be in the rotating frame and the velocity in the laboratory frame.
In the past, Bloch equation based MRI simulators were developed for variety of
research directions e.g. optimizing MR sequences, artifact detection, testing image reconstruction techniques, design of specialized RF pulses and educational purposes. Also,
multiple utilities of numerical simulations have been combined to produce a few general
purpose MRI simulators e.g., [5, 6, 7, 8, 9]. Accurate simulation of this initial-value
problem is still challenging for the following reasons: very tiny time steps, sufficiently
fine spatial resolution, and non-smooth data (e.g. gradient field G). To overcome this
difficulty, the execution speeds were improved using parallelization with message passing
interface (MPI), e.g. [8, 6], or GPU e.g. [9]. There are two prevailing techniques for the
development of the simulator (i) a semi-analytical technique based on operator splitting
as used, e.g., in [8, 9] (ii) coupled approach with fast ODE-solvers [6].
Similarly, there are many numerical studies on the influence of flow on MRI e.g., [10,
11, 4, 12]. Modified Bloch model was solved using multiple numerical strategies previously
e.g., Jou et al. [11], Lorthois et al. [4] first solved the flow field in a computational mesh
using finite volume method (FVM) softwares and then studied the effect of flow on
magnetization using finite difference method (FDM). Jurczuk et al. [12] solved 1.5 by
splitting the transport and the MR terms. In their work, the magnetization transport was
calculated using lattice-Boltzmann method (LBM) and the reaction part was calculated
using operator splitting techniques.
3

However, the previous studies lack a proper analysis of the Bloch model for flowing
objects. Sec. 2 is devoted to prove the well-posedness of the modified Bloch model. The
spatial semi-discretization with DG-methods is presented in Sec. 3 which, to the best of
our knowledge, is also an addition to the literature related to Bloch model. This section
was followed by consideration of the temporal discretization in Sec. 4. Verifications for
simple problems for static and flowing objects are presented in Sec. 5.
2. Well-posedeness of the Bloch Model
Let Ω ⊂ R3 be the flow domain with piecewise smooth Lipschitz boundary ΓSand outer
S
normal n. The incompressible flow field u(t, r) introduces the splitting Γ = Γ− Γ+ Γ0
where Γ− = {r ∈ Γ|u · n < 0}, Γ+ = {r ∈ Γ|u · n > 0} and Γ0 = {r ∈ Γ|u · n = 0}
represent the inflow boundary, outflow boundary and solid wall, respectively. We assume
that inflow and outflow are separated, i.e. dist(Γ− , Γ+ ) := min(P,Q)∈Γ− ×Γ+ |P − Q| > 0.
T
Problem (1.5) can be rewritten using Beff = B = Bx , By , Bz , the relaxation time

T
0
diagonal matrix D = diag( T12 , T12 , T11 ) and the constant source term f = 0, 0, M
T1
along with appropriate boundary and initial conditions as
∂M
+ (u · ∇)M + γB × M + DM = f ,
(t, r) ∈ [0, T ] × Ω,
∂t
M = MΓ , (t, r) ∈ [0, T ] × Γ− ,

(2.1a)
(2.1b)

0

M=M ,

(t, r) ∈ {0} × Ω.
(2.1c)
R
Consider the space H := [L2 (Ω)]3 with inner product (u, v)H := Ω u·v dr and norm
p
kvkH := (v, v). Moreover, define
X = {N ∈ H : (u · ∇)N ∈ H}

(2.2)

 12

kNkX := k(u · ∇)Nk2H + kNk2H .

(2.3)

with graph norm

Multiplying 2.1a by an arbitrary test function N ∈ X, integrating over Ω and imposing
(2.1b) weakly, we obtain
(∂t M, N)H + ((u · ∇)M, N)H + γ(B × M, N)H
(2.4)
Z
Z
+ (DM, N)H + (u · n)⊖ M · N ds = (f , N)H + (u · n)⊖ MΓ · N ds,
Γ

Γ

⊖

where w (r) :=

1
2 ( w(r)

− w(r)) and M × B = −B × M. Define bilinear and linear forms

a(t; M, N) := ((u · ∇)M, N)H + γ((B × M), N)H
Z
+ (DM, N)H + (u · n)⊖ M · N ds,
Z Γ
l(N) := (f , N)H + (u · n)⊖ MΓ · N ds,
Γ

4

(2.5a)
(2.5b)

Then we obtain the variational form of the Bloch problem: Find M : (0, T ] → X s.t.
(∂t M, N)H + a(t; M, N) = l(N),
0

M|t=0 = M .

∀ N ∈ X,

(2.6)
(2.7)

This is a Friedrichs system, see [13], Sec. 7. Unfortunately, the theory in [13, 14]
is not applicable since some coefficients are time-dependent. For 0 < ǫ ≪ 1, we set
Xǫ := [W 1,2 (Ω)]3 and consider an elliptic regularization of (2.6): Find Mǫ : (0, T ] → Xǫ
s.t.
(∂t Mǫ , N)H + aǫ (t; Mǫ , N) = l(N),
0

Mǫ |t=0 = M .

∀ N ∈ Xǫ ,

(2.8)
(2.9)

with
aǫ (t; M, N) := a(t; M, N) + ǫ(∇M, ∇N)H .

(2.10)

Please note that the variational formulation of the regularized Bloch problem incorporates
do-nothing boundary conditions ǫ∇Mǫ · n = 0 on Γ0 ∪ Γ+ .
The spaces Xǫ ⊆ H and the dual space Xǫ ∗ form an evolution triple (Xǫ , H, Xǫ ∗ ).
For p ≥ 1 and Banach space Y we denote by Lp (0, T ; Y) the Bochner space of vectorvalued functions v : (0, T ) → Y. We look for a solution Mǫ ∈ L∞ (0, T ; H) ∩ L2 (0, T ; Xǫ )
of problem (2.8)-(2.9).
Theorem 2.1 (Well-posedness). For all ǫ > 0, for given u ∈ [L∞ (0, T ; W 1,∞ (Ω)]3 with
div u = 0 and B ∈ [L∞ (0, T ; H)]3, there exists a unique solution Mǫ ∈ L∞ (0, T ; H) ∩
L2 (0, T ; Xǫ ) to Problems 2.8–2.9. For t ∈ (0, T ] and with σ := T11 , the following a-priori
estimate is valid


Z t
Z
1
1
2
2
σ(τ −t)
ǫ ∇Mǫ (τ ) H +
|(u · n)|(Mǫ · Mǫ )(s, τ ) ds dτ
Mǫ (t) H +
e
2
2 Γ
0
Z t
1
1
2
2
≤
Mǫ (0) H e−σt +
f (τ ) H eσ(τ −t) dτ.
(2.11)
2
2σ 0
Proof. In order to prove well-posedness, we apply the main existence theorem for evolution problems by J.L. Lions, see Theorem 6.6 in [13]. To this end, we define the norm
h
i1
2
2 2
kMkXǫ = ǫk∇MkL2 (Ω) +kMkX .
(2.12)

For the application of the theorem, the following conditions must be satisfied:

(P1) The time-dependent bilinear form t 7→ aǫ (t; M, N) is measurable for all M, N ∈ Xǫ
since the fields u and B are sufficiently smooth.
(P2) Application of the Cauchy-Schwarz and generalized Hölder inequalities show that
the bilinear form aǫ (t; ·, ·) is bounded for t ∈ [0, T ] and for all M, N ∈ Xǫ :
aǫ (t; M, N) ≤ǫk∇MkH k∇NkH + (u · ∇M, N)H + γ (B × M, N)H
Z
+ (DM, N)H + (u · n)⊖ M · N ds
Γ

≤ǫk∇MkH k∇NkH +ku · ∇MkH kNkH + γkBkL∞ kMkH kNkH
+kDkL∞ kMkH kNkH + ks kMkX kNkX
5

with kBk∞ := kBkL∞ (0,T ;[L∞ (Ω)]3 ) , kDkL∞ := max{ T11 , T12 } = T12 as T1 ≥ T2 and
where we used the trace inequality (see [14], Lemma 2.5)
Z
(u · n)⊖ M · N ds ≤ ks ku · ∇MkH ku · ∇NkH ≤ ks kMkX kNkX .
Γ

Now the norm definitions (2.3) and (2.12) imply boundedness:
aǫ (t; M, N) ≤(2 + γkBkL∞ +kDkL∞ + ks )kMkXǫ kNkXǫ .
(P3) Finally, bilinear form aǫ fulfills a coercivity condition. Integration by parts yields
Z
Z
2
1
1
aǫ (t; N, N) = ǫk∇Nk2H + D /2 N
+
(u · n)N · N ds + (u · n)⊖ N · N ds
2 Γ
H
Γ
Z
2
1
1
(u · n) N · N ds
= ǫk∇Nk2H + D /2 N
+
2
H
Z Γ
1
2
2
≥ ǫk∇NkH + σkNkH +
(u · n) N · N ds
(2.13)
2 Γ
with
1
1
1
1
D /2 := diag( √ , √ , √ ),
T2
T2
T1

σ := min(

1 1
1
, )=
(as T1 ≥ T2 ),
T1 T2
T1

(2.14)

as (B × N, N) = 0 and div u = 0 due to the incompressibility assumption.
Now we can apply the theorem by Lions giving existence and uniqueness of a generalized
solution Mǫ : [0, T ] 7→ X of (2.8)–(2.9).
It remains to prove the a-priori estimate. For all t ∈ [0, T ], we take N = Mǫ in
(2.8)–(2.9) to obtain via Cauchy-Schwarz and Young inequalities
1
σ
1 d
kMǫ k2H + aǫ (t; Mǫ , Mǫ ) = (f , Mǫ )H ≤
kf k2H + kMǫ k2H ,
2 dt
2σ
2
then via (P3)
1 d
σ
1
2
2
2
kMǫ kH + ǫk∇Mǫ kH + kMǫ kH +
2 dt
2
2

Z

Γ

(u · n) Mǫ · Mǫ ds ≤

1
2
kf kH .
2σ

Now, application of the Gronwall Lemma, see Lemma 6.9 in [13], implies


Z t
Z
1
1
2
kMǫ k2H +
eσ(τ −t) ǫ ∇Mǫ (τ ) H +
|(u · n)|(Mǫ · Mǫ )(s, τ ) ds dτ
2
2 Γ
0
Z t
1
1
2
2
≤
Mǫ (0) H e−σt +
f (τ ) H eσ(τ −t) dτ.
(2.15)
2
2σ 0
and we obtain estimate (2.11).

6

Remark. For the given case of constant source f , we obtain
Z
M0 2
M0 2
2
2
) dr = (
) |Ω| ,
f (t) H = kf kH = (
T
T1
1
Ω

|Ω| =

Z

dr

(2.16)

Ω

hence
Z

t

f (τ )

0

2 σ(τ −t)
e
dτ
H

=(

M0 2
1 − e−σt
) |Ω|
.
T1
σ

Finally, we can pass to the limit ǫ → +0, i.e. to the Bloch model.
Theorem 2.2. The Bloch model (2.6) admits a unique solution M ∈ L∞ (0, T ; H) ∩
L2 (0, T ; X). The kinetic energy of the magnetic field is bounded by:
Z

Z
1
1 t σ(τ −t)
2
e
M(t) H +
|(u · n)|(M · M))(s, τ ) ds dτ
2
2 0
Γ
Z t
1
1
2
2
M(0) H e−σt +
≤
f (τ ) H eσ(τ −t) dτ.
(2.17)
2
2σ 0
Proof. A careful inspection of the proof of Theorem 2.1 shows that the existence and
uniqueness result together with the a-priori estimate remain valid for ǫ → +0. As already
mentioned, the variational formulation of the regularized Bloch problem incorporates donothing boundary conditions ǫ∇Mǫ · n = 0 on Γ0 ∪ Γ+ . In the proof of the limit ǫ → +0
at Γ0 ∪ Γ+ , one can proceed as in Chapter V.1 of [15]. Note that here it is used that
inflow and outflow are separated.
Remark. The result of Theorem 2.1 and a-priori estimate 2.17 remain valid for the
special case u = 0, i.e. Bloch equations for spatially stationary objects.
3. Semi-discrete Equation
Here we consider the spatial discretization of problem (2.6) by dG-FEM [14].
3.1. Discontinuous Galerkin Formulation
Consider a non-overlapping decomposition Th := {Ωi }Ii=1 into convex simplicial subdomains Ωi , i = 1, 2, · · · , I as depicted in the left part of Fig. 2. We define the discontinuous finite element space
[Pk (Th )]d := {Nh ∈ H; Nh |Ωi ∈ [Pk (Ωi )]d

∀ Ωi , i = 1, 2, · · · , I}

(3.1)

T
where Pk denotes the set of polynomials of degree k ∈ N. Moreover,
let Xh = [Pk (Th )]d X.
T
For adjacent subdomains Ωi , Ωj with interface E = Γij = Ωi Ωj and unit normal vector
nij (directed from Ωi to Ωj ), as depicted in the right part of Fig. 2, we define the average
and jump of Nh ∈ X across Γij by
1
(Nh |Ωi (r) + Nh |Ωj (r)),
2
[Nh ]Γij (r) := Nh |Ωi (r) − Nh |Ωj (r).

hNh iΓij (r) :=

7

(3.2a)
(3.2b)

Figure 2: Left: 2D-simplicial mesh. Right: 1D-example of average and jump operators. Interface E
between i-th and j-th cell (j > i) is depicted with the used notation. Orientation of outward normal is
from lower to higher numbered cell.

Let Fhi be the set of all the interior interfaces E ⊆ Ω and define the upwind form

X Z 
1
− (u · nE )[M] · hNi + |u · nE | [M] · [N] ds.
Sh (t; M, N) :=
(3.3)
2
E
i
E∈Fh

Moreover, gradient jumps over interior faces are penalized via
Z
X
2
pǫ (M, N) := ǫ̃
hE
|u · nE |[∇M]E : [∇N]E ds,
E∈Fhi

E

ǫ̃ ≥ 0.

(3.4)

Setting
aupw
(t; M, N) := pǫ (M, N) + a(t; M, N) + Sh (t; M, N),
ǫ

(3.5)

the upwind dG-FEM reads: find Mh : (0, T ] 7→ Xh such that ∀ Nh ∈ Xh
(∂t Mh , Nh )H + aupw
(t; Mh , Nh ) = l(Nh ).
ǫ
3.2. Well-posedness of the Semi-discrete Problem
Let us define the norm |||Nh |||U via
Z
X
2
|||Nh |||U := ǫ̃
h2E
|u · nE | [∇Nh ]E
E∈Fhi

1
+
2

E

2
L2 (E)

Z
1 X
|u · n||Nh | ds +
|u · nE | [Nh ]2E ds.
2
Γ
E
i

Z

2

(3.6)

(3.7)

E∈Fh

Theorem 3.1. The semi-discrete problem 3.6 is well-posed and admits the a-priori
estimate
Z t
Z t
1
1
1
2
2
2
2
Mh (t) H +
eσ(τ −t) |||Mh (τ )|||U dτ ≤
Mh (0) H e−σt +
eσ(τ −t) f (τ ) H dτ.
2
2
2σ 0
0
(3.8)
8

Proof. The existence and uniqueness proof follows the lines of the proof of Theorem 2.1.
Similiar to the approach in Sec. 2, symmetric testing Nh = Mh provides
aupw
((t; Mh , Mh ) ≥ σkMh k2H + |||Mh |||2U .
ǫ

(3.9)

Then, Young inequality shows
1 d
Mh (t)
2 dt

2
H

+

σ
1
2
2
2
kMh kH + |||Mh |||U ≤
kf kH .
2
2σ

Now, similarly to Theorem 2.2, application of the Gronwall lemma yields the a-priori
estimate (3.8). This shows the well-posedness of the semi-discretized Bloch problem.
3.3. Semi-discrete Error Estimate
For the error of the spatial discretization we obtain the following result.
Theorem 3.2. The error of spatial discretization is given by:
Z
1
1 t σ(τ −t)
1
2
2
2
(M − Mh )(t)) H +
(M − Mh )(0) H e−σt
e
|||(M − Mh )(τ )|||U dτ ≤
2
2 0
2
Z t


2
+
eσ(τ −t) |||(M − πh M)(τ )|||U,♭ + δ (M − πh M)(τ ) H + pǫ (M, M)(τ ) dτ (3.10)
0

with the L2 -orthogonal projection πh M of M onto Xh and σ = T11 . The norm |||·|||U,♭
and constant δ will be defined within the proof. For a sufficiently smooth solution M ∈
L∞ (0, T ; [W k+1,2 (Ω)]3 ), the right-hand side term in (3.10) is of order O h2k+1 with
h = maxi diam(Ωi ).
Proof. The error equation for the error M − Mh is given by
(∂t (M − Mh ), Nh )H + aupw
(t; M − Mh , Nh ) = pǫ (M, Nh ) ∀ Nh ∈ Xh , t ∈ (0, T ] a.e..
ǫ
Let πh N be the L2 -orthogonal projection of N onto Xh , i.e., (N− πh N, w)H = 0 ∀ w ∈
Xh . Then we split the error as
M − Mh = (M − πh M) + (πh M − Mh ) ≡ Ih + Eh
and reformulate the error equation with Nh = Eh as
(∂t Eh , Eh )H + aupw
(t; Eh , Eh ) = −aupw
(t; Ih , Eh ) + pǫ (M, Eh ),
ǫ
ǫ

t ∈ (0, T ] a.e. (3.11)

where we used that (∂t Ih , Eh )H = 0 due to L2 -orthogonality. The left-hand side can be
bounded from below as
(t; Eh , Eh ) ≥
(∂t Eh , Eh )H + aupw
ǫ

1 d
2
2
2
kEh kH + σkEh kH + |||Eh |||U .
2 dt

(3.12)

For an estimate of the right-hand side we use pǫ (M, Eh ) ≤ pǫ (M, M)|||Eh |||U and


2
−aupw
(t;
I
,
E
)
≤
|||I
|||
+
δkI
k
(3.13)
h
h
h
h
ǫ
U,♭
H |||Eh |||U
9

with
2
2
|||Ih |||U,♭ := ζ|||Ih |||U +

X

T ∈Th

2

kukL∞ (∂T ) kIh kL2 (∂T )

(3.14)

and δ := γkBkL∞ +kDkL∞ resp. ζ := max{1;kukL∞ (0,T ;W 1,∞ (Ω)) }. Estimate (3.13) relies
on a generalization of Lemma 2.30 in [14] for scalar advection-reaction problems to the
vector-valued case. The estimate provides a careful bound of the uwpind-discretized
convective term and heavily exploits the boundedness of L2 -orthogonality of subscales.
Combining (3.12), (3.13) and using Young inequality, we obtain
1 d
1
2
2
2
2
kEh kH + |||Eh |||U ≤ |||Ih |||U,♭ + δkIh kH + pǫ (M, M).
2 dt
2

(3.15)

Integration of (3.15) and the triangle inequality imply the quasi-optimal error estimate
(3.10). For a sufficiently smooth solution M ∈ L∞ (0, T ; [W k+1,2 (Ω)]3 ), the penalty
 term
pǫ (M, M) vanishes. Finally, interpolation results imply the error order O hk+1/2 .

Remark. The case k = 0 covers the finite volume method (FVM).
4. Temporal Discretization
Starting point is the spatially discretized problem: find Mh : (0, T ] 7→ Xh such that
for all Nh ∈ Xh
(∂t Mh (t), Nh )H + aupw
(t; Mh (t), Nh ) = l(Nh ),
ǫ

Mh (0) = Mh0 .

(4.1)

A major problem stems from the multiscale character of (4.1) as the scale of magnetization is much faster than that of advection. Another difficulty is the restricted smoothness
T
of the data in time, in particular of field G = Gx , Gy , Gz , see Fig. 3 (for the FLASH
sequence [16, 17]).
According to the required high resolution in time, an explicit time stepping is chosen.
We considered two variants (i) a fully coupled approach and (ii) an operator splitting
approach. The efficiency of the numerical simulation can be strongly improved using
GPU computing. This will be exemplarily shown in Sec. 5.
4.1. Fully Coupled Approach
Following Sec. (3.1) in [14], we apply a low-order explicit Runge-Kutta scheme to (4.1).
: X + Xh 7→ Xh via (Aupw
(t; v, w).
(t)v, w)H := aupw
Define a discrete operator Aupw
ǫ
ǫ
ǫ
Similarly, let L be a functional on Xh with L = l(w). Note that L is constant in time in
this application.
Let 0 = t0 < t1 < t2 < · · · < tN = T be the set of discrete times with time steps
τn := tn+1 − tn , n = 0, 1, · · · , N − 1. Moreover, we denote Mnh = Mh (tn ) etc.
A two-stage RK-scheme is a good compromise between temporal accuracy and the
restricted data smoothness in time. Following Subsec. 3.1.3 in [14] we select scheme
n
upw
Mn,1
Mnh + τn L
h = Mh − τn Aǫ
1
1
1
upw
Mn,1
Mn+1
= (Mnh + Mn,1
h ) − τn Aǫ
h + τn L.
h
2
2
2
10

(4.2a)
(4.2b)

Precise statements of the stability and convergence of 4.2a can be found in case of smooth
data in Subsec. 3.1.6 of [14]. In particular, a time step restriction on τn comes from a
CFL condition for the advective term. As in this application the time step restriction
on τn comes from magnetization, the mentioned CFL condition is always valid in our
calculations.
We will not repeat the details, e.g. of the stability and convergence
 in
 RK2-analysis
2
k+1/2
with
[14]. It provides in case of smooth data in time, an error of order O τn + h
polynomial degree k of spatial discretization (see Theorem 3.10). Such error estimate in
time is not valid in this application, since the data B(t, r) are only in [C 0,1 [0, T ]]3 (for
the FLASH-sequence studied in Sec. 5) or even only in [L∞ (0, T )]3 for the example in
Subsec. 5.2. Alternatively, a standard embedded RK-scheme of type RK3(2) or even of
higher order like RK5(4) is chosen for appropriate time step selection, see Sec. 5. For
details of the embedded methods, we refer, e.g., to Sec. 5 of [18].

4.2. Operator Splitting
The above mentioned multiscale character of problem (4.1) suggests an operator splitting approach as suggested e.g.in Sec. IV.1.5. of [19] for large advection-reaction problems, e.g. in air-pollution simulations. Writing the semi-discrete problem (4.1) with
Fadv (t, M) := −∇ · (uM),

Fmag := γB × M + DM − f ,

(4.3)

formally as ODE-system
dMh (t)
= Fadv (t, Mh (t)) + Fmag (t, Mh (t)),
dt

(4.4)

the simplest sequential operator splitting on tn ≤ t ≤ tn+1 gives
dM∗h (t)
= Fadv (t, M∗h (t)),
dt
dM∗∗
h (t)
= Fmag (t, M∗∗
h (t)),
dt

M∗h (t) = Mh (tn ),

(4.5a)

∗ n
M∗∗
h (t) = Mh (t ).

(4.5b)

An inspection of 4.5 shows for the splitting error at t = tn
ǫs =

 
1 2 ∂ Fadv
∂ Fmag
τn [
Fmag −
Fadv ] + O τn3 .
2
∂M
∂M

It turns out that the commutation error [Fadv , Fmag ] =
be written as
[Fadv , Fmag ] = (∇ · u)[Fmag (t, Mh ) −

(4.6)

∂ Fmag
∂ Fadv
Fmag −
Fadv can
∂Mh
∂Mh

∂ Fmag
(t, Mh )Mh (t)] + (u · ∇r )Fmag (t, Mh ).
∂t
(4.7)

It vanishes if either Fmag is independent of r = (x, y, z)T and div u = 0 or Fmag is
independent of r and linear in Mh , see [19], sec IV.1.5. This is unfortunately not the
case in this application. As a remedy, symmetric
 or Strang-Marchuk splitting can be
applied which reduces the splitting error to O τn2 , see [19].
11

Echo
RF
d

c
d

c
a
b
TE
TR
Figure 3: Generic spoiled gradient echo sequence with radial trajectory. Gradients: (a) slice selection
(b) rewinder (c) prephasing (d) readout. The dashed radial line (spoke) in the k-space corresponds to
the current repetition.

5. Numerical Experiments
The simulation method was tested for both static objects, i.e. with u ≡ 0, as well as
flow experiments, i.e. with u 6≡ 0.
5.1. Experimental Validation for Static Objects (u ≡ 0)

We start with a test of the simulator for the static case, i.e. for u = 0. This means
that mainly the temporal discretization is considered. The experiments here and later
on (with exception of Subsec. 5.2) were performed using a randomly spoiled [20] FLASH
sequence [16, 17] with a radial trajectory as shown in Fig. 3. The left part depicts the
time diagram of gradients and the right part shows the k-space sampling trajectory where
Z t
Z t
γ
γ
Gx (τ ) dτ,
ky :=
Gy (τ ) dτ
kx :=
2π 0
2π 0
q
with Gx = Gmax cos(θ), Gy = Gmax sin(θ), Gmax = G2x + G2y and angle θ of the spoke
with x-axis. Putting the expression of k in the demodulated signal equation, we can
obtain the Fourier pair relation between the acquired electric signal and the transverse
magnetization [1]. Each radial spoke in k-space, which represents a projection of the
object, is acquired with a repetition of pulse sequence where transverse gradients are
changed according to orientation θ of the spoke in k-space. Certain number of spokes in
k-space are required to reconstruct an image frame.
Carr showed in [21] that under constant flip angle α, and gradient moment and
constant TR the magnetizations reach a state of dynamic equilibrium after several repetitions. For clinical imaging, the acquisition starts only after the magnetizations reach
dynamic equilibrium after several preparatory TR repetition. In order to test the simulations, time-series data of images from the beginning of the preparatory phase to dynamic
equilibrium were obtained from experiments and compared with equivalent simulation
results.
12

To this end, an experiment was performed with Table 1: Doped water tubes with their rea phantom as can be seen in the left part of Fig. 4) laxation times
containing multiple compartments of doped water
Tube T1 [ms] T2 [ms]
tubes with known T1 and T2 as listed in Table 1.
3
296
113
The relevant pulse sequence parameters in Fig. 3
4
463
53
are: TR/TE = 2.18/1.28 ms, flip angle α =8◦ , and
7
604
95
number of spokes per image frame= 27. 100 image
10
745
157
frames from the beginning to the dynamic equilib14
1034
167
rium of the experiment were used for comparison
16
1276
204
with simulation results.
water
2700
2100
Before embarking on the validation of our code
with experimental results, it was tested with a Bloch equation simulator written in MATLAB by Sun [22] where the computationally expensive parts are written in C and connected via mex interface [22]. The code is an extended and corrected version of the code
by Hargreaves [7]. The numerical results from our code agree really well with [22]. Baseline simulations were performed on a system with Supermicro SuperServer 4027GR-TR
system with Ubuntu 14.04, 2x Intel Xeon Ivy Bridge E5-2650 main processors using an
in-house C++ code and the implementation by Sun in Matlab R2016A. Though computationally RK5(4) solver is 4 times more expensive, the advantage of using RK5(4) is an
easy higher order extension for flowing cases with any method of our choice for spatial
discretization and we can avoid splitting error using a fully coupled approach.
Prior studies [23] have shown that, to model Bloch phenomena accurately, simulations
must be carried out dividing each voxel with a number of cells where each of them
represent one or an ensemble of isochromats; an isochromat is a microscopic group of
spins which resonate at the same frequency. Shkarin et. al. showed in that [23] increasing
the number of cells per voxel reduce the error of the solution.
For every simulation with relaxation time of specific liquids, the simulated data were
recorded at TE since all isochromats are rephased at TE under the condition of approximately complete spoiling of residual transverse magnetization. The transverse magnetizations are integrated over all the isochromats and the magnitude of the summed up
transverse magnetizations are averaged further for the number of spokes per image frame.
In principle, this averaged integrated transverse magnetization should be equivalent to
the averaged magnitude of image over a region in a specific tube.
Simulations were performed over a domain of 4.8 × 4.8 × 18.0 mm3 (corresponding
to 3 × 3 pixels in the xy-plane and three-times the nominal slice thickness=6 mm in zdirection for each tube) divided into 27 × 27 × 45 cells. Each cell is assumed to consist
of one isochromat as shown in Subsec. 3.4.3. in [24] that the number of isochromats
in each cell do not effect accuracy of the simulation. The initial condition is chosen
as M0 = (0, 0, 1)T . The embedded RK5(4) scheme is applied for time discretization.
The time series of these two equivalent quantities are plotted as a function of image
frame number in Fig. 4 (right) for four randomly chosen tubes after normalizing the
experimental and simulated data by their respective magnitude in dynamic equilibrium
of the brightest tube (Tube 3).
The excellent agreement between simulation and experiment hints at the possible use
of the simulator for quantitative estimations in MRI, e.g. the relaxation times T1 , T2
or concentration of contrast agents required for certain signal enhancement (illustrated
with an example in [24], Sec. 5.4).
13

1.8
Experiment tube3
Experiment tube4
Experiment tube10
Experiment water tube
Simulation with relaxation times of the liquid in tube 3
Simulation with relaxation times of the liquid in tube 4
Simulation with relaxation times of the liquid in tube 10
Simulation with relaxation times of water

1.6
1.4

14

W
7
3

|M xy | [a.u]

10
16

1.2
1
0.8

4
0.6
0.4
0.2

10

20

30

40

50

60

70

80

90

100

Frame Number

Figure 4: Left: Image of the static phantom with resolution of 160 × 160 pixels. Right: Comparison of
simulation with the image for four different liquids.

Also, the time evolution of magnetizations for different isochromats are independent
of each other and suitable for Graphical Processing Unit (GPU) computing. A speed-up
of 82 was achieved with GPU parallelization on the previously mentioned system together
with a NVIDIA GTX Titan Black (Kepler GK110) GPU as illustrated in Subsec. 3.3.4
of [24].
5.2. A Basically One-dimensional Test Case for Flowing Spins
The simulator is tested further for flowing spins, i.e. u 6≡ 0. As first case, we consider
a basically one-dimensional test case in [10] where they studied for z ∈ (− L2 , L2 ) the effect
of an RF pulse on the magnetization for different through-plane velocities, i.e. component
uz in z-direction, using the leap-frog finite difference scheme.
As shown in Fig. 5, a Blackman-windowed sinc
pulse with an amplitude of 0.1750 G, flip angle of
90◦ and duration of 2.6794 ms, slice selection gradient Gz = 1.0 G cm−1 , and a nominal slice thickness of 7 mm were used for the simulation in [10].
Through-plane velocities were in the range of 0 to
200 cm s−1 . The simulations were performed over a
length of L = 20 mm and L = 30 mm, respectively,
in
the slice direction z, divided into 800 cells of size
◦
Figure 5: A 90 slice-selective pulse was
used for the studying the flow-effects. The 0.025 mm and 0.0375 mm for the range of velocities
arrows indicate the time when the data 0 to 80 cm s−1 and 80 to 200 cm s−1 respectively.
was recorded.
The time duration of simulations was divided into
4500 equidistant time steps of τn = 8.9313 × 10−4 ms. The magnetizations were calculated at the end of the rewinder gradient marked by an arrow in Fig. 5
In this work, the simulations were performed with the FEM-package COMSOL Multiphysics using the dG-FEM with quadratic elements (k = 2), see Sec. 3. The domain
of L=30 mm was divided into 200 equidistant cells for the velocity range 0 to 200 cm s−1 .
The fully coupled semi-discretized equations were further solved using the RK5(4) scheme.
14

For relative and absolute tolerance of 1.0 × 10−6 and 1.0 × 10−8 respectively, the adaptive time stepping resulted in 340 to 528 time steps in the specified velocity range. The
variability of time step sizes due to the embedded RK is depicted in Fig. 6 (Left) and the
figure also clearly shows sudden reduction in the time step size due to time adaptivity
to cater for jump in Gz . Also, the number of time steps got reduced by an order of
magnitude due to time adaptation of the embedded RK scheme.
In order to estimate grid convergence, simulations were performed for 16, 32, 64, 128
and 256 equidistant grid points in z-direction with uz = 80 cm s−1 . The solution with the
finest grid was chosen as the reference for the estimation of errors for My . Fig. 6 (Right)
shows the error order and the results show faster convergence rate than the predicted
theoretical estimate.

10

8

10

7

10

-2

10 -4

10 6

10

O(h 2 )
O(h 2.5 )
O(h 3 )
H 1 Error
L 2 Error

10 -6

Error

Reciprocal of step size

10 0

10

-8

10

-10

10

-12

10

-14

5

10 4

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time [ms]

-3.7

-3.6

-3.5

-3.4

-3.3

-3.2

-3.1

-3

-2.9

-2.8

-2.7

log10(h)

Figure 6: Left: Reciprocal of variable time step size as a function of time duration for uz = 200 cm s−1 .
Right: Numerical errors estimated with different grid sizes comparing with the reference solution obtained dividing the length in z-direction into 256 grid points and with uz = 80 cm s−1 . Different error
estimates are plotted in logarithmic scale as a function of logarithm of grid sizes.

The penalty diffusion term pǫ with constant ǫ̃ = 5.0 × 10−4 was used to remove
unphysical oscillations appearing in the neighborhood of inflow boundary for Mx and
My . Such oscillations for My are shown in Fig. 7. Similar oscillations were observed for
Mx but not for Mz .
Finally, the results are compared with [10] and they show extremely good agreement.
The effect of velocity uz in the range 0 to 200 cm s−1 on magnetization is shown in Fig. 8.
5.3. Comparison with Experiments for Through-plane Flow
The simulator was evaluated further against a laminar flow experiment in a circular
tube. In the experiment, the flow pump was operated at different voltages to produce
velocities such that the flow profile could be expected to be laminar i.e. Reynolds number
Re ≤ 2300 as listed in Table 2.
The flow velocities were estimated pixelwise using phase contrast MRI (PC MRI) [25].
At each listed operating voltage in Table 2, the mean through-plane flow velocity uz was
calculated over a pixel. The measured velocities show an unsteady pattern with a mean
and a standard deviation as listed in the second and the third column of Table 2.

15

Figure 7: Simulation for uz = 160 cm s−1 . Left: with ǫ = 0 results in unphysical oscillations in the
neighborhood of inflow boundary for My . Right: Unphysical oscillations vanish after the introduction
of a small artificial diffusion term with ǫ̃ = 5.0 × 10−4 .
Table 2: Mean and standard deviation velocities and Re based on the mean velocity for different operating
voltages of the flow pump at temperature 16 ◦C (kinematic viscosity ν = 1.1092 × 10−2 cm2 /s).

Voltage [V]
6
5
4
3

Mean Velocity [mm s−1 ]
49.19
38.71
28.84
18.52

Standard Deviation [mm s−1 ]
2.26
1.97
1.47
1.04

Re
2217
1744
1300
834

At each specified flow velocity uz , measurements were performed with the FLASH
pulse sequence parameters TR/TE = 1.96/1.22 ms, flip angle = 8◦ , and number of spokes
per image frame = 17. The resolution of one pixel in the xy-plane is 1.6 × 1.6 mm2
whereas the nominal slice thickness in z-direction is 6 mm. The time series of averaged
magnitude is recorded over the same region of interest as the region of flow velocity
measurement. 60 frames from the beginning of experiments towards dynamic equilibrium
were used for the comparison with simulation.
The mean velocity u = (0, 0, uz )T was taken as input velocity in the simulation. A
computational domain of 1.6 × 1.6 × 18.0 mm3 , divided into 6 × 6 × 27 identical cells,
was chosen for the simulation. The length of domain in the through-plane flow direction
was estimated from previous simulation results [24] such that MΓ− = (0, 0, 1)T , i.e. the
Dirichlet inflow boundary condition, could be satisfied. The initial condition was chosen
as in the static case. The spatial discretization was performed using dG-FEM with
quadratic elements (k = 2) whereas in time the fully coupled approach with the RK5(4)scheme for time discretization was applied. Again the penalty diffusion term pǫ with
constant ǫ = 5.0 × 10−4 was used.
For comparison, the experimental and simulated data are normalized properly first.
Due to unsteadiness of the flow, the experimental data was normalized by the average
magnitude of last 20 frames from the measurement with the operating voltage of 6 V in
Table 2. The simulated data was normalized similarly. After that, they are plotted in
the right part of Fig. 9.
16

0.8
80

120

0.4
200

Mx 0
0

-0.4

-0.8
160
-1.2
1

120

0.5
0

My

0
160

-0.5

80

200

-1
1.2

1
0

0.8

80

M z0.6

120
160

0.4

200

0.2
0
-0.2

-10

-5

0

5

10

Z [mm]

Figure 8: Simulated magnetization distributions of Mx , My , Mz for through-plane velocity uz along the
positive z-axis in the range 0 to 200 cm s−1 using dG-method (right) are compared with the results in
[10] (Left). The magnetizations were recorded at the end of rewinder gradient as marked by the arrow
Fig. 5. The length in the slice direction is in the range −15 to 15 mm.

In spite of the unsteadiness in the flow which may be attributed to the lack of sufficient
entry length, as discussed in [24], the simulation and the experiment show a reasonable
agreement.
Like the static phantom, GPU computing can be used for the flowing fluid as well to
achieve a significant speed up as explained in Subsec. 4.4.3 of [24].

17

1.3
1.2

u z [mm s -1 ]
49.19
2.26
38.71
1.97
28.84
1.47
18.52
1.04
49.19
38.71
28.84
18.52

|M xy | [a.u]

1.1
1
0.9
0.8
0.7
0.6

5

10

15

20

25

30

35

40

45

50

55

60

Frame Number

Figure 9: Left: MR image of the flow tube. Right: Experiments with different velocities are compared
with the simulations.

5.4. Comparison with Pulsatile Flow Experiments
The simulation method was evaluated further with a pulsatile flow laboratory experiment. Like in the previous case, the pulsatile velocity profile was estimated using PC
MRI. A through-plane profile was fitted as a function of time using Matlab curve-fitting
toolbox as shown in the left part of Fig. 10. The fitted pulsatile flow profile was used as
input velocity in the simulation. The effect of pulsatile flow on the evolution of magnitude from the image was studied. An experiment was performed using the pulse sequence
identical to the previous experiment.
10

2.5

9

2

7

|M xy | [a.u]

u z [mm s -1 ]

8

6
5

1.5

4

Experiment
Simulation

3
2

2

3

4

5

6

7

8

1

9

Time [s]

2

3

4

5

6

7

8

9

Time [s]

Figure 10: Left: Fitted through-plane velocity profile from PC MRI data. Right: Comparison between
experimental and simulation results for studying the effect of pulsatile flow.

The simulation was performed with identical domain, spatial and temporal discretization, and penalty diffusion as used in the previous experiment. Experimental data was
normalized with the magnitude data in dynamic equilibrium from the spatially stationary
tube under identical experimental conditions. The simulated data was normalized by the
averaged magnitude of integrated transverse magnetization of static water in dynamic
equilibrium. The experiment and the simulation are compared in a state of dynamic
equilibrium in Fig. 10 (right).
18

Fig. 10 shows that, although the periodicity in the experimental and the simulated
result agree well, the amplitude of experimental results deviates from the simulation. The
deviation could be due to the assumption of a flow profile which depends only on time.
The assumption implies that the fluid needs to move in bulk, i.e., the tube must respond
simultaneously to the changing pressure at all positions at every specific point of time in
the direction of through-plane flow such that through-plane velocities at every position
in the longitudinal direction are same, which is artificial and unphysical. Nevertheless,
it is a starting point to study the effect of more realistic pulsatile flows.
6. Summary. Outlook
In this note, we proved the well-posedness of the Bloch model under the action of
an incompressible flow field. Then we applied the discontinuous Galerkin method to the
spatial semi-discretization of the Bloch model and proved well-posedness and error estimates. The multiscale character of the problem basically stems from the high frequency
time evolution of the magnetization part. An explicit Runge-Kutta method together with
time step adaption is applied for the temporal discretization. Alternatively, an operator
splitting between advection and magnetization can be applied. The computation can be
strongly accelerated via GPU computing.
Magnetic resonance imaging is nowadays a very rapid process which can be done
in real-time [26]. Nevertheless, there are still unsolved problems such as a quantitative
understanding of the mechanisms that lead to MRI signal alterations (i.e., both enhancement and loss) when imaging flowing spins (e.g., in vessels or the heart) or other dynamic
processes. Here, numerical simulation with the proposed direct solver for the Bloch model
can help in a better understanding of such dynamic processes as shown for some basic
MRI experiments. Apart from flow velocities and volumes, there is an increasing demand
in MRI for quantitative information such as relaxation time constants. In future, access
to both high-contrast imaging and quantitative parametric mapping by MRI is expected
to facilitate and contribute to computer-aided diagnostic strategies.
Acknowledgement
The first author thanks Prof. Dr. J. Frahm of the Biomedizinische NMR Forschungs
GmbH (BiomedNMR) for suggesting the topic of MRI and many helpful discussions and
Dr. D. Voit (BiomedNMR) for stimulating discussions on various occasions and numerous
guidance for the experiments. We thank P. Schroeder for his helpful remarks. Moreover,
we would also like to thank two anonymous referees for their helpful comments.
References
[1] D. G. Nishimura, Principles of magnetic resonance imaging, Stanford University, 1996.
[2] F. Bloch, Nuclear induction, Physical Review 70 (7-8) (1946) 460–474.
[3] E. O. Stejskal, J. E. Tanner, Spin diffusion measurements: spin echoes in the presence of a timedependant field gradient, The Journal of chemical physics 42 (1) (1965) 288–292.
[4] S. Lorthois, J. Stroud-Rossman, S. Berger, L. D. Jou, D. Saloner, Numerical simulation of magnetic resonance angiographies of an anatomically realistic stenotic carotid bifurcation, Annals of
Biomedical Engineering 33 (3) (2005) 270–283.

19

[5] J. Bittoun, J. Taquin, M. Sauzade, A computer algorithm for the simulation of any Nuclear Magnetic
Resonance (NMR) imaging method, Magnetic Resonance Imaging 2 (2) (1984) 113–120.
[6] T. Stöcker, K. Vahedipour, D. Pflugfelder, N. J. Shah, High-performance computing MRI simulations, Magnetic Resonance in Medicine 64 (1) (2010) 186–193. doi:10.1002/mrm.22406.
[7] B. Hargreaves, Bloch equation simulator, http://mrsrl.stanford.edu/~ brian/blochsim/, published January 13, 2004.
[8] H. Benoit-Cattin, G. Collewet, B. Belaroussi, H. Saint-Jalmes, C. Odet, The SIMRI project: A
versatile and interactive MRI simulator, Journal of Magnetic Resonance 173 (1) (2005) 97–115.
[9] C. G. Xanthis, I. E. Venetis, A. V. Chalkias, A. H. Aletras, MRISIMUL: A GPU-based parallel
approach to MRI simulations, IEEE Transactions on Medical Imaging 33 (3) (2014) 607–617.
[10] C. Yuan, G. T. Gullberg, D. L. Parker, The solution of bloch equations for flowing spins during a
selective pulse using a finite difference method, Medical physics 14 (6) (1987) 914–921.
[11] L. D. Jou, R. van Tyen, S. A. Berger, D. Saloner, Calculation of the magnetization distribution for
fluid flow in curved vessels, Magnetic Resonance in Medicine 35 (4) (1996) 577–584.
[12] K. Jurczuk, M. Kretowski, J. J. Bellanger, P. A. Eliat, H. Saint-Jalmes, J. Bézy-Wendling, Computational modeling of MR flow imaging by the lattice Boltzmann method and Bloch equation,
Magnetic Resonance Imaging 31 (7) (2013) 1163–1173.
[13] A. Ern, J.-L. Guermond, Theory and practice of finite elements, Vol. 159, Springer Science &
Business Media, 2004.
[14] D. A. Di Pietro, A. Ern, Mathematical aspects of discontinuous Galerkin methods, Vol. 69, Springer
Science & Business Media, 2011.
[15] J. L. Lions, Perturbations singulières dans les problèmes aux limites et en contrôle optimal, Vol.
323, Springer, 2006.
[16] J. Frahm, A. Haase, D. Matthaei, W. Hänicke, K.-D. Merboldt, Verfahren und Einrichtung zur
schnellen Akquisition von Spinresonanzdaten für eine ortsaufgelöste Untersuchung eines Objekts
(1985).
[17] A. Haase, J. Frahm, D. Matthaei, W. Hanicke, K.-D. Merboldt, FLASH imaging. Rapid NMR
imaging using low flip-angle pulses, Journal of Magnetic Resonance 67 (2) (1986) 258–266.
[18] P. Deuflhard, F. Bornemann, Scientific computing with ordinary differential equations, Vol. 42,
Springer Science & Business Media, 2012.
[19] W. Hundsdorfer, J. G. Verwer, Numerical Solution of Time-Dependent Advection-DiffusionReaction Equations, Springer Science & Business Media, 2003.
[20] V. Roeloffs, D. Voit, J. Frahm, Spoiling without additional gradients: Radial flash mri with
randomized radiofrequency phases, Magnetic Resonance in Medicine 75 (5) (2016) 2094–2099.
doi:10.1002/mrm.25809.
[21] H. Y. Carr, Steady-state free precession in nuclear magnetic resonance, Physical Review 112 (5)
(1958) 1693–1701.
[22] H. Sun, Bloch equation simulator, https://github.com/ismrmrd/ismrmrd-paper/blob/master/code/extern/irt/mri-rf/sun-bl
published Jul 22,2012.
[23] P. Shkarin, R. G. S. Spencer, Time domain simulation of Fourier imaging by summation of isochromats, International journal of imaging systems and technology 8 (1997) 419–426.
[24] A. Hazra, Numerical Simulation of Bloch Equations for Dynamic Magnetic Resonance Imaging,
Ph.D. thesis, Institut für Numerische und Angewandte Mathematik (2016).
[25] A. Joseph, J. T. Kowallick, K. D. Merboldt, D. Voit, S. Schaetz, S. Zhang, J. M. Sohns, J. Lotz,
J. Frahm, Real-time flow MRI of the aorta at a resolution of 40 msec, Journal of Magnetic Resonance
Imaging 40 (1) (2014) 206–213.
[26] M. Uecker, S. Zhang, D. Voit, A. Karaus, K.-D. Merboldt, J. Frahm, Real-time mri at a resolution
of 20 ms, NMR in Biomedicine 23 (8) (2010) 986–994.

20

