Doron Levy

We adopt a non-oscillatory  central  scheme, first presented in the context of Hyperbolic conservation laws in [nessyahu-tadmor] followed by [jiang-tadmor], to the framework of the incompressible Euler equations in their vorticity formulation. The embedded duality in these equations, enables us to toggle between their two equivalent representations - the conservative Hyperbolic-like form vs. the convective form. We are therefore able to apply local methods, to problems with a global nature. This results in a new stable and convergent method which enjoys high-resolution without the formation of spurious oscillations. These desirable properties are clearly visible in the numerical simulations we present. 

The integration of semi-discrete approximations for time-dependent problems is encountered in a variety of applications. The Runge-Kutta (RK) methods are widely used to integrate the ODE systems which arise in this context, resulting in large ODE systems called methods of lines. These methods of lines are governed by possibly ill-conditioned systems with a growing dimension: consequently, the naive spectral stability analysis based on scalar eigenvalues arguments may be misleading. Instead, we present here a stability analysis of RK methods for well-posed semi-discrete approximations, based on a general energy method. We review the stability question for such RK approximations, and highlight its intricate dependence on the growing dimension of the problem. In particular, we prove the strong-stability of general fully-discrete RK methods governed by coercive approximations. We conclude with two non-trivial examples which demonstrate the versatility of our approach in the context of general systems of convection-diffusion equations with variable coefficients. A straightforward implementation of our results verify the strong-stability of RK methods for local finite-difference schemes as well as global spectral approximations. Since our approach is based on the energy method (which is carried in the physical space), and since it avoids the von-Neumann analysis (which is carried in the dual Fourier space), we are able to easily adapt additional extensions due to non-periodic boundary conditions, general geometries, etc. 

We present a general procedure to convert schemes which are based on staggered spatial grids into non-staggered schemes. This procedure is then utilized to construct a new family of non-staggered, central schemes for hyperbolic conservation laws, by converting the family of staggered central schemes recently introduced in [nessyahu-tadmor], [liu-tadmor], [jiang-tadmor]. These new non-staggered central schemes retain the desirable properties of simplicity and high-resolution, and in particular, they yield Riemann-solver-free recipes which avoid dimensional splitting. Most importantly, the new central schemes avoid staggered grids and hence are simpler to implement in frameworks which involve complex geometries and boundary conditions. 

We analyze self-focusing and singularity formation in the complex Ginzburg-Landau equation (CGL) in the regime where it is close to the critical nonlinear Schrodinger equation. Using modulation theory [Fibich and Papanicolaou, Phys. Lett. A 239:167--173, 1998], we derive a reduced system of ordinary differential equations that describes self-focusing in CGL. Analysis of the reduced system shows that in the physical regime of the parameters there is no blowup in CGL. Rather, the solution focuses once and then defocuses. The validity of the analysis is verified by comparison of numerical solutions of CGL with those of the reduced system. 

We present a new central scheme for approximating solutions of two-dimensional systems of hyperbolic conservation laws. This method is based on a modification of the staggered grid proposed in [jiang-tadmor] which prevents the crossings of discontinuities in the normal direction, while retaining the simplicity of the central framework. Our method satisfies a local maximum principle which is based on a more compact stencil. Unlike the previous method, it enables a natural extension to adaptive methods on structured grids. 

The KdV equation can be considered as a special case of the general equation $ u_{t} + f(u)_{x} - \delta g(u_{xx})_x = 0, \delta > 0$, where $f$ is non-linear and $g$ is linear, namely $f(u)=u^2/2$ and $g(v)=v$. As the parameter $\delta$ tends to $0$, the dispersive behavior of the KdV equation has been throughly investigated (see, e.g., [whitham], [lax-levermore], [drazin] and the references therein). We show, through numerical evidence that a completely different, dissipative behavior occurs when $g$ is non-linear, namely when $g$ is an even concave function such as $g(v)=-|v|$ or $g(v)=-v^2$. In particular, our numerical results hint that as $\delta$ tends to zero, the solutions strongly converge to the unique entropy solution of the formal limit equation, in total contrast with the solutions of the KdV equation.  

We consider a family of high-order, weighted essentially non-oscillatory central schemes (CWENO) for approximating solutions of one-dimensional hyperbolic systems of conservation laws. We are interested in the behavior of the total variation (TV) of the approximate solution obtained with these methods. Our numerical results suggest that even though CWENO methods are not total variation diminishing (TVD), they do have bounded total variation (TVB). Moreover, the TV of the approximate solution seems to never increase above the theoretical value, and it approaches it as the mesh is refined. These results are hopefully a first step in the quest for proving the convergence of such high-order methods. 

We present a new third-order central scheme for approximating solutions of systems of conservation laws in one and two space dimensions. In the spirit of Godunov-type schemes, our method is based on reconstructing a piecewise-polynomial interpolant from cell-averages which is then advanced exactly in time. In the reconstruction step, we introduce a new third-order,  compact, CWENO reconstruction, which is written as a convex combination of interpolants based on different stencils. The heart of the matter is that one of these interpolants is taken as an arbitrary quadratic polynomial and the weights of the convex combination are set as to obtain third-order accuracy in smooth regions. The embedded mechanism in the WENO-like schemes guarantees that in regions with discontinuities or large gradients, there is an automatic switch to a one-sided second-order reconstruction, which prevents the creation of spurious oscillations. In the one-dimensional case, our new third order scheme is based on an extremely compact  four  stencil. Analogous compactness is retained in more space dimensions. The accuracy, robustness and high-resolution properties of our scheme are demonstrated in a variety of one and two dimensional problems. 

We introduce a new dispersion-velocity particle method for approximating solutions of linear and nonlinear dispersive equations. This is the first time in which particle methods are being used for solving such equations. Our method is based on an extension of the diffusion-velocity method of Degond and Mustieles [SIAM J. Sci. Stat. Comp., 11 no. 2, (1990), pp.293--310] to the dispersive framework. The main analytical result we provide is the short time existence and uniqueness of a solution to the resulting dispersion-velocity transport equation. We numerically test our new method for a variety of linear and nonlinear problems. In particular we are interested in nonlinear equations which generate structures that have non-smooth fronts. Our simulations show that this particle method is capable of capturing the nonlinear regime of a compacton-compacton type interaction. 

We present the first fourth-order, central scheme for two-dimensional hyperbolic systems of conservation laws. Our new method is based on a Central Weighted Non-Oscillatory (CWENO) approach. The heart of our method is the reconstruction step, in which a genuinely two-dimensional interpolant is reconstructed from cell-averages by taking a convex combination of building blocks in the form of bi-quadratic polynomials. Similarly to other central schemes, our new method enjoys the simplicity of the black-box approach. All that is required in order to solve a problem is to supply the flux function and an estimate on the speed of propagation. The high-resolution properties of the scheme as well as its resistance to mesh orientation, and the effectiveness of the component-wise approach, are demonstrated in a variety of numerical examples. 

We present new third- and fifth-order Godunov-type central schemes for approximating solutions of the Hamilton-Jacobi (HJ) equation in an arbitrary number of space dimensions.

We present new third- and fifth-order Godunov-type central schemes for approximating solutions of the Hamilton-Jacobi (HJ) equation in an arbitrary number of space dimensions. These are the first central schemes for approximating solutions of the HJ equations with an order of accuracy that is greater than two. In two space dimensions we present two versions for the third-order scheme: one scheme that is based on a genuinely two-dimensional Central WENO reconstruction, and another scheme that is based on a simpler dimension-by-dimension reconstruction. The simpler dimension-by-dimension variant is then extended to a multi-dimensional fifth-order scheme. Our numerical examples in one, two and three space dimensions verify the expected order of accuracy of the schemes. 

We develop local discontinuous Galerkin (DG) methods for solving nonlinear dispersive partial differential equations that have compactly supported traveling waves solutions, the so-called "compactons". The schemes we present extend the previous works of Yan and Shu on approximating solutions for linear dispersive equations and for certain KdV-type equations. We present two classes of DG methods for approximating solutions of such PDEs. First, we generate nonlinearly-stable numerical schemes with a stability condition that is induced from a conservation law of the PDE. An alternative approach is based on constructing linearly-stable schemes, i.e. schemes that are linearly stable to small perturbations. The numerical simulations we present verify the desired properties of the methods including their expected order of accuracy. In particular, we demonstrate the potential advantages of using DG methods over pseudo-spectral methods in situations where discontinuous fronts and rapid oscillations co-exist in a solution. 

Recently, a wavelet-based method was introduced for the systematic derivation of subgrid scale models in the numerical solution of partial differential equations. Starting from a discretization of the multiscale differential operator, the discrete operator is represented in a wavelet space and projected onto a coarser subspace. The coarse (homogenized) operator is then replaced by a sparse approximation to increase the efficiency of the resulting algorithm. In this work we show how to improve the efficiency of this numerical homogenization method by choosing a different compact representation of the homogenized operator. In two dimensions our approach for obtaining a sparse representation is significantly simpler than the alternative sparse representations. $L^{\infty}$ error estimates are derived for a sample elliptic problem. An additional improvement we propose is a natural fine-scales correction that can be implemented in the final homogenization step. This modification of the scheme improves the resolution of the approximation without any significant increase in the computational cost. We apply our method to a variety of test problems including one- and two-dimensional elliptic models as well as wave propagation problems in materials with subgrid inhomogeneities. 

We introduce a new family of Godunov-type semi-discrete central schemes for multidimensional Hamilton-Jacobi equations. These schemes are a less dissipative generalization of the central-upwind schemes that have been recently proposed in [Kurganov, Noelle and Petrova, SIAM J. Sci. Comput., 23 (2001), pp. 707-740]. We provide the details of the new family of methods in one, two and three space dimensions, and then verify their expected low-dissipative property in a variety of examples. 

In inverse planning, the likelihood for the points in a target or sensitive structure to meet their dosimetric goals is generally heterogeneous and represents the a priori knowledge of the system once the patient and beam configuration are chosen. Because of this intrinsic heterogeneity, in some extreme cases, a region in a target may never meet the prescribed dose without seriously deteriorating the doses in other areas. Conversely, the prescription in a region may be easily met without violating the tolerance of any sensitive structure. In this work, we introduce the concept of dosimetric capability to quantify the a priori information and develop a strategy to integrate the data into the inverse planning process. An iterative algorithm is implemented to numerically compute the capability distribution on a case specific basis. A method of incorporating the capability data into inverse planning is developed by heuristically modulating the importance of the individual voxels according to the a priori capability distribution. The formalism is applied to a few specific examples to illustrate the technical details of the new inverse planning technique. Our study indicates that the dosimetric capability is a useful concept to better understand the complex inverse planning problem and an effective use of the information allows us to construct a clinically more meaningful objective function to improve IMRT dose optimization techniques. 

We present two new model equations for the unidirectional propagation of long waves in dispersive media for the specific purpose of modeling water waves. The derivation of the new equations uses a Pade (2,2) approximation of the phase velocity that arises in the linear water wave theory. Unlike the Korteweg-deVries (KdV) equation and similarly to the Benjamin-Bona-Mahony (BBM) equation, our models have a bounded dispersion relation. At the same time, the equations we propose provide the best approximation of the phase velocity for small wave numbers that can be obtained with third-order equations. We note that the new model equations can be transformed into previously studied models, such as the BBM and the Burgers-Poisson equations. It is therefore straightforward to establish the existence and uniqueness of solutions to the new equations. We also show that the distance between the solutions of one of the new equations, the KdV equation, and the BBM equation, is of the small order that is formally neglected by all models. 

We derive a second-order, semi-discrete central-upwind scheme for the incompressible two-dimensional Euler equations in the vorticity formulation. The reconstructed velocity field preserves an exact discrete incompressibility relation. We state a local maximum principle for a fully-discrete version of the scheme and prove it using a convexity argument. We then show how similar convexity arguments can be used to prove that the scheme maps certain Orlicz spaces into themselves. The consequences of this result on the convergence of the scheme are discussed. Numerical simulations support the expected properties of the scheme. 

A multiscale image registration technique is presented for the registration of medical images that contain significant levels of noise. An overview of the medical image registration problem is presented, and various registration techniques are discussed. Experiments using mean squares, normalized correlation, and mutual information optimal linear registration are presented that determine the noise levels at which registration using these techniques fails. Further experiments in which classical denoising algorithms are applied prior to registration are presented, and it is shown that registration fails in this case for significantly high levels of noise, as well. The hierarchical multiscale image decomposition of E. Tadmor, S. Nezzar, and L. Vese \cite{tadmor} is presented, and accurate registration of noisy images is achieved by obtaining a hierarchical multiscale decomposition of the images and registering the resulting components. This approach enables successful registration of images that contain noise levels well beyond the level at which ordinary optimal linear registration fails. Image registration experiments demonstrate the accuracy and efficiency of the multiscale registration technique, and for all noise levels, the multiscale technique is as accurate as or more accurate than ordinary registration techniques.

We incorporate new high-order WENO-type reconstructions into Godunov-type central schemes for Hamilton-Jacobi equations. We study schemes that are obtained by combining the Kurganov-Noelle-Petrova flux with the Weighted Power ENO and the Mapped WENO reconstructions. We also derive new variants of these reconstructions by composing the Weighted Power ENO and the Mapped WENO reconstructions with each other. While all schemes are, formally, fifth-order accurate, we show that the quality of the approximation does depend on the particular reconstruction that is being used. In certain cases, it is shown that the approximate solution may not converge to the viscosity solution at all. 

This paper focuses on the characterization of delay effects on the asymptotic stability of some continuous-time delay systems encountered in modeling the post-transplantation dynamics of the immune response to chronic myelogenous leukemia. More explicitly, we shall discuss the stability of the crossing boundaries of the corresponding linearized models in the delay-parameter space. Weak, and strong cell interactions are discussed, and analytic characterizations are proposed. An illustrative example completes the presentation. 

In this paper we introduce new models for describing the motion of phototaxis, i.e. bacteria that move towards light.   Following experimental observations, the first model describes the locations of bacteria, an internal property of the bacteria that is related to the group dynamics, and the interaction between the bacteria and the medium in which it resides.   The second model is a new multi-particle system following the same quantities as in the first model.  The main theorem shows how to obtain a new system of PDEs, which we refer to as the phototaxis system as the limit dynamics of the multi-particle system.  Numerical simulations are provided.

An image registration technique is presented for the registration of medical images using a hybrid combination of coarse-scale landmark and B-splines deformable registration techniques.  The technique is particularly effective for registration problems in which the images to be registered contain large localized deformations.  A brief overview of landmark and deformable registration techniques is presented.  The hierarchical multiscale image decomposition of E. Tadmor, S. Nezzar, and L. Vese, "A multiscale image representation using hierarchical $(BV,L^2)$ decompositions", Multiscale Modeling and Simulations, vol. 2, no.4, pp. 554-579, 2004, is reviewed, and an image registration algorithm is developed based on combining the multiscale decomposition with landmark and deformable techniques. Successful registration of medical images is achieved by first obtaining a hierarchical multiscale decomposition of the images and then using landmark-based registration to register the resulting coarse scales.  Corresponding bony structure landmarks are easily identified in the coarse scales, which contain only the large shapes and main features of the image.  This registration is then fine-tuned by using the resulting transformation as the starting point to deformably register the original images with one another using an iterated multiscale B-splines deformable registration technique.  The accuracy and efficiency of the hybrid technique is demonstrated with several image registration case studies in two and three dimensions.  Additionally, the hybrid technique is shown to be very robust with respect to the location of landmarks and presence of noise. 

This paper focuses on the characterization of delay effects on the asymptotic stability of some continuous-time delay systems encountered in modeling the post-transplantation dynamics of the immune response to chronic myelogenous leukemia. Such models include multiple delays in some large range, from one minute to several days. The main objective of the paper is to discuss the stability of the crossing boundaries of the corresponding linearized models in the delay-parameter space by taking into account the interactions between small and large delays. Weak, and strong cell interactions are discussed, and analytic characterizations are proposed. An illustrative example together with related discussions complete the presentation.

This work presents a hierarchy of mathematical models for describing the motion of phototaxis, i.e., bacteria that move towards light.  Based on experimental observations, we conjecture that the motion of the colony towards light depends on certain group dynamics.  This group dynamics is assumed to be encoded as an individual property of each bacterium, which we refer to as 'excitation'.  The excitation of each individual bacterium changes based on the excitation of the neighboring bacteria. Under these assumptions, we derive a stochastic model for describing the evolution in time of the location of bacteria, the excitation of individual bacteria, and a surface memory effect.   A discretization of this model results in an interacting stochastic many-particle system.  The third, and last model is a system of partial differential equations that is obtained as the continuum limit of the stochastic particle system.  The main theoretical results establish the validity of the new system of PDEs as the limit dynamics of the multi-particle system.

We derive a model for describing the dynamics of imatinib-treated chronic myelogenous lekemia (CML).  This model is a continuous extension of the agent-based CML model of Roeder et al.  [Nature Medicine, 12(10), 1181-1184, 2006] and of its recent formulation as a system of difference equations [Kim, Lee, Levy, Bull. Math. Bio., 70, 2008 728-744].  The new model is formulated as a system of partial differential equations that describe various stages of differentiation and maturation of normal hematopoietic cells and of leukemic cells.  An imatinib treatment is also incorporated into the model.  The simulations of the new PDE model are shown to qualitatively agree with the results that were obtained with the discrete-time (difference equations and agent-based) models.  At the same time, for a quantitative agreement, it is necessary to adjust the values of certain parameters, such as the rates of imatinib-induced inhibition and degradation. 

Adaptive radiation therapy (ART) is the incorporation of daily images in the radiotherapy treatment process so that the treatment plan can be evaluated and modified to maximize the amount of radiation does to the tumor while minimizing the amount of radiation delivered to healthy tissue.  Registration of planning images with daily images is thus an important component of ART.  In this article, the authors report their research on multiscale registration of planning computed tomography (CT) images with daily cone beam CT (CBCT) images.  The multiscale algorithm is based on the hierarchical multiscale image decomposition of E. Tadmor, S. Nezzar, and L. Vese [Multiscale Model. Simul. 2(4), pp. 554-579 (2004)].  Registration is achieved by decomposing the images to be registered into a series of scales using the (BV,L^2) decomposition and initially registering the coarsest scales of the image using a landmark-based registration algorithm.  The resulting transformation is then used as a starting point to deformably register the next coarse scales with one another.  This procedure is iterated at each stage using the transformation computed by the previous scale registration as the starting point for the current registration.  The authors present the results of studies of rectum, head-neck, and prostate CT-CBCT registration, and validate their registration method quantitatively using synthetic results in which the exact transformations our known, and qualitatively using clinical deformations in which the exact results are not known.

We study a stochastic Cucker-Smale flocking system in which particles interact with the environment through white noises. We provide definition of flocking to the stochastic system, and show that when the communication rate is constant, the system exhibits a flocking behavior independent of the initial configurations. For the case of a radially symmetric communication rate with a positive lower bound, we show that the relative fluctuations of the particle velocity around the mean velocity have a uniformly bounded variance in time. We conclude with numerical simulations that validate our analytical results.